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ABSTRACT 


M 


A  nonlinear  controller  is  designed  for  a  full  size 
submarine  using  the  LQG/LTR  procedure.  Linear  models  of  the 
submarine  are  developed  at  four  different  speeds  and  then 
analyzed  using  the  method  of  modal  analysis.  The  linear 
models  are  then  augmented  with  integral  control  and  a  Kalman 
Filter  transfer  function  is  designed  using  some  tools  for 
loop  shaping.  The  Loop  Transfer  Recovery  technique  is  then 
applied  to  recover  the  Kalman  Filter  loop  shapes.  A 
polynomial  data  fit  is  performed  on  the  resulting 
compensators  to  produce  a  nonlinear  controller.  Both  the 
linear  and  the  nonlinear  controllers  are  extensively  tested 
using  a  full  nonlinear  model  of  the  submarine, 
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Chapter  1 


INTRODUCTION  AND  8UHMARY 


1.1  Background 

The  technology  for  building  multivariable  control  systems 
has  advanced  in  the  last  several  years  and  there  is  a  very 
definite  need  for  complex  design  examples  to  help  engineers 
understand  and  further  develop  this  powerful  design 
methodology. 

Multi-input,  multi-output  (MIMO)  control  system  design  is 
much  more  complex  than  either  classical  control  system 
design  or  single-input,  single-output  (SISO)  control  system 
design.  The  engineer  must  use  and  manipulate  the  singular 
values  of  the  loop  transfer  function  matrix  of  the  plant  to 
determine  the  performance  of  the  controller  and  this  concept 
of  singular  values  as  a  measure  is  at  present  relatively 
unknown  to  most  control  engineers.  Therefore,  it  is 
necessary  to  apply  the  MIMO  methodology  to  practical 
examples  in  order  to  display  the  power  of  this  technique  and 


B 


also  to  point  out  any  shortcomings. 


To  date,  non-trivial  examples  of  MIMO  design  using  the 
Linear  Quadratic  Gaussian  procedure  with  Loop  Transfer 
Recovery  (LQG/LTR)  are  few  in  number  (see  Cl 3  C 21  for  recent 
examples).  This  is  in  part  due  to  the  LQG/LTR  technique 
having  been  only  recently  developed  C33  C43  C5D ,  and  partly 
due  to  the  significant  effort  required  to  develop  and 
linearize  a  MIMO  model,  design  the  controller,  and  test  the 
resulting  design. 


Previous  controller  designs  for  submarines  have  used  the 
SISO  design  technique  or  classical  design  techniques  CAJ 
Z71,  This  thesis  is  believed  to  be  the  first  (unclassified) 
example  of  a  MIMO  design  on  a  full-size  submarine.  The 
LQG/LTR  design  methodology  is  used  because  it  incorporates 
the  dynamics  of  the  plant  (the  submarine)  into  the 
compensator  and  thereby  provides  a  vehicle  for  capturing  the 
dynamics  of  the  submarine  at  various  speeds.  This 
information  can  then  be  utilized  to  produce  a  nonlinear 
dynamic  compensator  that  varies  as  some  sensed  parameter 
changes  (gain  scheduling),  providing  adaptive-like  nonlinear 
control  of  the  submarine  over  a  variety  of  speed  conditions 
and  simultaneous  maneuvers.  In  fact,  the  linearized  models 
of  the  submarine  were  selected  so  as  to  emphasize  the 
cross-coupling  of  the  longitudinal  and  lateral  dynamics. 
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1.2  Contributions  of  the  Thesis 


The  main  contribution  o-f  this  thesis  is  to  illustrate  the 
multi  variable  LQG/LTR  -feedback  control  system  design 
methodology  -for  a  maneuvering  submarine,  and,  further,  to 
demonstrate  the  utilization  of  gain  scheduling  to  construct 


a  nonlinear  controller.  The  reader  shoul d  not  interpret  the 
resulting  control ler  design  to  be  a  final  product ,  and 
further  research  is  needed  to  use  such  controllers  for  all 
diverse  submarine  maneuvers. 


The  first  stage  of  the  design  process  is  to  achieve  as 
complete  an  understanding  of  the  submarine  model  as  is 
possible.  To  accomplish  this,  the  submarine  model  is 
linearized  about  an  unconventional  operating  trajectory  so 
as  to  capture  longi tudinal /lateral  interactions.  The 
resulting  linear  representation  is  decomposed  into  its  modal 
representation.  The  modal  matrices  are  then  graphically 
displayed  to  visualize  the  dynamic  behavior  of  the  submarine 
and  further  manipulated  to  determine  whether  or  not  the 
required  observabi 1 i ty/control labi 1 i ty  requirements  of  the 
LQG/LTR  design  process  are  met  under  various  choices  of 
control  and  output  variables. 
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The  controller  design  is  based  on  the  LQ6/LTR  methodology 
formulated  by  Doyle  and  Stein  C33  C43,  and  refined  by  Athens 


The  design  effort  concentrates  on  the  frequency  domain 
properties  of  the  state  space  representations  of  the 
submarine  model.  Methods  of  shaping  the  singular  values  of 
the  open  loop  plant  transfer  matrix  are  developed  that 
guarantee  the  identical  behavior  of  the  singular  values  at 
both  the  high  and  low  ends  of  the  frequency  spectrum. 
Scaling  of  the  output  variables  is  utilized  to  produce  the 
desired  loop  shapes  at  frequencies  near  crossover. 


The  state  variable  used  for  gain  scheduling  was  taken  to 
be  forward  velocity  in  the  body  reference  coordinate 
system.  Although  gain  scheduling  may  also  be  desired  based 
on  roll  angle  (again  in  the  body  reference  coordinate 
system),  this  would  necessarily  be  in  addition  to  the  gain 
scheduling  employed  using  forward  velocity.  Time  domain 
plots  of  nonlinear  simulations  of  the  compensated  system 
designed  at  a  particular  speed  are  presented  to  show  the 
validity  of  that  choice. 


Gain  scheduling  of  the  compensator  designs  is  accomplished 
by  performing  a  least  squares  fit  on  the  individual 
compensator  parameters.  An  algorithm  for  producing  the 
coefficients  resulting  from  a  quadratic  polynomial  least 


squares  -fit  to  matrices  is  presented. 

Finally,  time  domain  plots  of  the  gain  scheduled  LQG/LTR 
compensator  being  applied  to  the  nonlinear  submarine  model 
are  presented  as  a  means  of  illustrating  the  resulting 
nonlinear  design. 


1.3  Outline  of  the  Thesis 


Chapter  2  describes  the  process  used  to  model  the 
submarine  and  the  methods  used  to  produce  a  linear  state 
space  representation.  The  computer  implementation  of  the 
submarine  nonlinear  and  linear  models  is  briefly  described. 
The  chapter  ends  with  a  discussion  of  the  reasoning  used  to 
select  the  output  variables  and  presents  the  performance 
requirements  imposed  on  the  controller  designs. 

In  Chapter  3,  the  eigenstructure  of  the  linear  model  is 
studied  using  modal  decomposition.  Pole-zero  information 
and  singular  value  plots  are  also  utilized  to  further 
display  the  structure  of  the  open  loop  dynamics. 

Chapter  4  contains  the  linear  portion  of  the  design, 
following  the  LQG/LTR  methodology.  Although  the  reader  is 
assumed  to  be  familiar  with  the  work  of  Doyle  and  Stein  C33, 
the  notation  being  used  is  briefly  summarized.  Singular 
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value  plots  of  the  open  loop  submarine  model  plus 

compensator  are  presented  for  each  design. 

The  gain  scheduling  method  for  the  compensator  parameters 
is  presented  in  Chapter  5.  Results  of  various  nonlinear 
simulations  are  included  to  display  the  performance, 
focusing  upon  simultaneous  evasive  maneuvers. 

Chapter  6  contains  a  summary,  conclusions,  and  some 
proposed  directions  for  further  research. 
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Chapter  2 


MATHEMATICAL  MODEL  OF  THE  SUBMARINE 


2. 1  Introduction 


The  submarine  model  used  in  this  thesis  is  that 
implemented  at  the  Charles  Stark  Draper  Laboratory  (CSDL) 
both  as  a  real-time  simulation  -facility  and  an  analytical 
model  generating  facility.  Details  of  the  latter  are  given 
in  [93 . 

This  chapter  discusses  the  equations  used  to  model  a 
submarine  and  presents  a  brief  history  of  the  implementation 
of  the  submarine  model  at  the  CSDL. 

The  software  in  its  present  form  Cl  13  3  provides  the  control 
system  designer  with  several  useful  tools.  The  manner  in 
which  these  tools  are  used  to  produce  and  validate  the 
required  linear  models  for  a  gain  scheduling  LQB/LTR 
controller  is  detailed,  as  well  as  the  reasoning  used  to 
select  the  output  variables. 
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The  chapter  ends  with  a  section  on  performance 
requirements  to  be  imposed  on  the  controller  design,  thus 
setting  the  stage  for  analysis  of  the  model  in  chapter  3  and 
finally  the  actual  design  of  the  controller  in  chapter  4. 


2.2  Modeling  of  the  Submarine 


The  generic  model  from  which  most  submarine  models  are 
derived  is  that  given  in  the  original  Naval  Ship  Research 
and  Development  Center  (NSRDC)  Report  2510  CUD  entitled 
"Standard  Equations  of  Motion  for  Submarine  Simulation." 
The  model  used  in  this  thesis  is  believed  to  be  an 


improvement  of  the  'standard  model '  in  that  it  includes  the 
crossflow  drag  and  vortex  contribution  terms  for  the  five 
degrees  of  freedom  (the  straight  ahead  x-axis  force  is  not 


included) . 


To  establish  the  notation  for  describing  the  submarine 
motion,  a  brief  description  of  the  nonlinear  equations  of 
motion  is  presented. 


Two  reference  frames  are  used  in  deriving  the  equations  of 


motion  for  a 


-  the  body  reference  frame  and  the 


earth  reference  frame. 


former  is  used  because 


hydrodynamic  forces  and  inertias  are  more  readily  computed 
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in  the  submarine  reference  frame.  For  purposes  of  control 
and  stability  determination t  me  are  also  concerned  with  the 
motion  of  the  submarine  in  the  earth  reference  frame.  The 
relationship  of  the  motion  of  the  submarine  between  these 
two  reference  frames  is  described  by  Euler's  angles  C121  V, 
0,  and  ♦* 

V  (Yaw)  -  rotation  about  the  z  axis 
0  (Pitch)  -  rotation  about  the  y  axis 
♦  (Roll)  -  rotation  about  the  x  axis 

where  x,  y,  and  z  represent  the  body  reference  frame  as  it 
changes  according  to  the  indicated  rotations  and  in  the 
order  given.  The  orientation  of  the  submarine  in  its  own 
reference  frame  is  depicted  in  figure  2.1. 


Figure  2.1*-  Submarine  Body  Reference  Frame  Axes 


The  state  vector  for  the  submarine  must  therefore  include 
the  six  degrees  of  freedom  in  its  own  reference  frame,  the 
three  Euler  angles,  and  any  desired  position  variables  to 

locate  the  submarine  with  respect  to  some  earth  frame 

reference  point.  For  the  purposes  of  this  thesis,  only  the 
depth  position  variable  is  used  and  the  state  vector  for  the 
submarine  model  contains  the  states  summarized  in  table  2.1. 

Table  2.1.  Submarine  State  Variables 
u  =  Xj  (t)  -  forward  velocity  in  feet /sec 

v  =  x_(t)  -  lateral  velocity  in  feet/sec 

w  =  x^(t)  “  vertical  velocity  in  feet/sec 
p  *  x^(t)  --  roll  rate  in  radians/sec 
q  *  _  pitch  rate  in  radians/sec 

r  =  x&<t)  “  yaw  rate  in  radians/sec 
♦  *  x7<t>  -  roll  angle  in  radians 
9  *  Xg(t)  -  pitch  angle  in  radians 
V  =  x9<t>  -  yaw  angle  in  radians 
z  *=  Xjg(t)  -  depth  in  feet  (positive  downward) 

Notet  u,  v,  w,  ♦  ,  and  9  are  in  the  body 
reference  frame.  All  others  are  in  the  inertial 
reference  frame. 

It  should  be  noted  that  the  Euler  angles  only 
approximately  describe  the  submarine  roll,  pitch,  and  yaw 
angles  in  the  earth  reference  frame,  with  the  approximation 
becoming  more  accurate  as  the  magnitude  of  the  Euler  angles 


approaches  zero.  This  is  due  to  the  fact  that,  in 
transforming  the  coordinates  from  earth  reference  to  body 
reference,  the  Euler  angle  rotations  are  applied 
sequentially,  in  the  specific  order  V,  0,  and  finally  ♦. 

Using  the  state  variables  of  table  2.1,  the  nonlinear 
equations  of  motion  are  derived  using  force  and  moment 
balances  in  the  submarine  rotating  coordinate  system.  The 
forces  and  moments  due  to  the  submarine  motion  must  be  equal 
to  the  total  forces  and  moments  exerted  by: 

1.  Hydrodynamic  pressures 

2.  Control  forces  and  moments 

3.  Propulsion  or  tow  forces 

4.  Any  other  forces  and  moments 

The  reader  is  referred  to  the  NSROC  Report  2510  1121  for 
details  of  the  derivation  of  the  nonlinear  equations  of 
motion  and  a  description  of  the  hydrodynamic  coefficients 
describing  the  submarine  geometry  and  control  surfaces. 

The  submarine  model  used  in  this  thesis  has  three  control 
variables  -  rudders,  stern  planes,  and  fairwater  planes  (see 
figure  2.2).  The  propeller  is  constrained  to  turn  at  a 
constant  rpm  to  reflect  current  operating  procedures.  The 
control  surfaces  on  the  submarine  are  locked  together  in  the 


controller  design,  the  same  set  of  nonlinear  equations  was 
later  implemented  on  the  IBM  time  sharing  computer  at  CSDL. 
Considerable  enhancement  was  added  to  the  program  during  the 
transfer  to  aid  the  control  system  designer  in  his  task.  At 
present,  capabilities  of  the  software  include! 

1.  Addition  of  a  user-f riendly  executive  routine  to  allow 
modification  of  parameters  and  selection  of  options 
for  a  particular  simulation  run.  This  routine  then 
submits  the  tailored  program  for  background 
execution. 

2.  The  option  of  calculating  A  and  B  matrices  that 
describe  the  linearization  of  the  submarine  about  a 
particular  nominal  point,  in  the  form 

*<t>  =  A  x(t)  +  B  u  ( t ) . 

3.  The  options  of  setting  control  surfaces  to  fixed 
values,  varying  the  control  surfaces  over  time 
according  to  values  in  a  data  file,  calculation  of  the 
controls  using  full  state  feedback,  or  calculation  of 
the  controls  using  a  LQ6/LTR  derived  compensation. 

4.  Selection  of  computer  print-out  or  plots  (or  both)  of 
the  state  variables  over  time  of  either  the  nonlinear 
model  or  the  linearized  model. 

5.  The  capability  of  searching  for  a  local  equilibrium 


point  for  the  nonlinear  model  that  is  close  to  a 
desired  nominal  point. 

To  closely  simulate  typical  operating  profiles  of  an 
actual  submarine,  the  propeller  rpm  is  kept  constant  during 
a  particular  run. 

2.4  Generation  of  the  Linear  Models 

To  perform  gain  scheduling  of  compensator  designs  over  a 
large  speed  range  requires  that  several  linear  models  of  the 
submarine  be  developed.  Extensive  analysis  of  the  submarine 
open  loop  dynamics,  however,  requires  only  a  single  model 
that  captures  most  of  the  dynamic  modes  present  in  the 
submarine.  As  we  remarked  before,  this  linearized  model  was 
obtained  to  capture  significant  longitudinal -lateral 
interactions  in  simultaneous  maneuvers. 

The  submarine  modeled  in  this  thesis  (figure  2.2)  is 
approximately  400  feet  long,  has  the  conventional  stern 
conf iguration  (stern  control  surfaces  are  in  the  shape  of  a 
cross),  and  planes  mounted  on  the  sail  (fairwater  planes). 
The  lack  of  differential  control  for  the  control  surfaces 
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reason,  the  initial  linear  model  should  be  such  that  roll 
does  not  dominate  the  submarine  behavior  nor  is  it  absent 
either.  The  nominal  point  used  for  analysis  uses  a  forward 
velocity  of  20  knots,  5  degrees  of  rudder,  and  one  degree  on 
each  of  the  planes.  This  places  the  submarine  in  a  diving 
turn  with  the  trajectory  that  of  a  helix.  Other  state 
variables  are  as  determined  by  an  equilibrium  point  search 
with  these  initial  conditions.  Analysis  of  this  model  is 
presented  in  the  next  chapter. 
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Figure  2.2.  Submarine  Control  Surfaces 


The  control  variables  used  in  the  controller  designs  are 
ds,  db,  and  dr  (refer  to  table  2.2),  the  angular  deflections 
of  the  three  control  surfaces.  Deflection  limits,  as 
imposed  by  the  software,  are  40  degrees  for  the  stern  planes 
and  rudder,  and  20  degrees  for  the  sail  planes. 


Linearization  of  the  nonlinear  equations  of  motion  is 


performed  in  a  straightf orward  manner.  A  nominal  point  is 
determined  by  integrating  the  nonlinear  equations  of  motion 
using  a  selected  set  of  initial  condi tions,  and  an 
equilibrium  point  is  found  that  corresponds  to  minimum 
accelerations  for  all  state  variables  very  near  the  steady 
state  values  determined  from  the  integration  of  the 
equations  of  motion.  For  the  nonlinear  equations  of  motion, 
there  are  11  state  variables  since  the  control  surface 
deflections  are  also  independent  variables  <¥  and  z  have  no 
effect  on  the  integration).  Thus  an  equilibrium  point  is 
more  precisely  defined  as  a  local  minimum  in  11-space.  The 
values  of  the  state  variables  at  the  equilibrium  point  then 
describe  a  nominal  point,  about  which  high  order  terms  can 
be  neglected.  A  set  of  linear  differential  equations  can 
then  be  produced,  and  the  A  and  B  matrices  calculated,  to 
provide  a  state  space  description  of  the  submarine. 

For  each  nominal  point  thus  determined,  the  resulting 
linear  model  must  be  validated  by  perturbing  the  nominal 
point  to  form  a  set  of  initial  conditions,  and  then 
comparing  the  results  of  integrating  both  the  nonlinear  and 
the  linear  equations  of  motion.  Provided  the  perturbations 
are  not  too  large,  the  nonlinear  model  will  always  return  to 
the  equilibrium  point  values,  while  the  linear  model  will 
never  reach  steady  state  due  to  the  non-zero  forcing  term 


imposed  by  the  control  surfaces.  The  comparison  of  the  two 
time  histories,  therefore,  is  limited  to  checking  initial 
derivatives,  apparent  natural  frequencies,  and  the  damping 
rate. 


Four  linear  models  were  developed,  corresponding  to 
forward  velocities  of  5,  10,  20,  and  25  knots.  Control 
surface  deflections  were  maintained  at  the  values  previously 
mentioned.  The  models  are  designated  S5R5,  S10R5,  S20R5, 
and  S25R5,  reflecting  the  speed  and  rudder  deflection. 
Figures  2.3  and  2.4  show  results  of  the  comparisons  for 
selected  state  variables  for  the  5  and  20  knot  models.  The 
perturbations  applied  to  the  nominal  point  were  identical  in 
all  cases  (see  table  2.3).  The  comparisons  show  excellent 
correlation  between  the  1 i near  and  nonlinear  models  and 
serve  to  validate  the  linear  models.  In  particular,  the 
initial  derivatives,  the  apparent  natural  frequency  of 
response,  and  the  damping  factors  match  almost  exactly  for 
the  linear  and  nonlinear  models.  Notice  that  the  pitch  rate 
q(t)  (figure  2.3  (c)  and  (d)  and  2.4  (c)  and  (d)>  contains 
two  oscillatory  modes,  reflecting  the  cross  coupling  between 
the  rudders  and  stern  planes  when  the  submarine  has  a 
non-zero  roll  angle.  The  roll  angle  for  the  5  knot  model  is 
-0.6  degrees  and  for  the  20  knot  model,  -8  degrees.  The  A 
and  B  matrices  for  the  four  models  are  presented  in  appendix 
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(e)  Yaw  rate  r(t)  for  the  linear  model 


Table  2.3 


Perturbations  Applied  to  Noainal  Points 


State  Variable 

u 

v 

M 

P 

q 

r 

♦ 

6 


Perturbation 

+2.0  ft/sec 
+0.2  ft/sec 

-  0.02  -ft/sec 

-  0.00003  rad/sec 
+  0.00003  rad/sec 

-  0.003  rad/sec 

-  0. S  degrees 

-  1.0  degrees 


2.5  Selection  of  the  Output  Variables 

Selection  of  the  output  variables  requires  both  a  careful 
study  of  the  A  and  B  matrices  and  determination  of  the 
intent  of  the  controller  design.  Three  control  variables 
are  available  so  that  three  output  variables  can  be  chosen. 
An  autopilot  could  be  designed f  using  the  position  variables 
V  and  or  a  rate  controller  could  be  designed,  using  the 
rate  variables  u,  v,  w,  p,  q,  or  r.  The  attitude  variables 
♦  and  O  could  be  included  in  either  design,  depending  on  the 
importance  of  these  variables  to  the  control  system 
designer. 

The  philosophy  taken  in  this  thesis  is  that  a  rate 
controller  be  designed  to  control  the  submarine  during 
various  maneuvers.  As  previously  mentioned,  the  submarine 
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model  was  linearized  about  a  nonconvent ional  operating  point 
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to  ensure  that  the  cross  coupling  of  the  control  surfaces 
Mould  be  captured  in  the  A  and  B  matrices.  The  simplifying 
assumption  is  made  that  the  rate  commands  are  presented  in 
the  form  of  commands  for  m  and  r.  If  this  assumption  Mere 
not  made,  the  coordinate  transf or mat ion  effect  of  the  Euler 
angles  Mould  require  a  C  matrix  that  Mas  dependent  on  the 
real  time  values  of  0,  and  y,  making  the  C  matrix  state 
dependent.  The  third  output  variable  is  taken  to  be  0,  the 
pitch  angle,  permitting  the  controller  to  accomplish  depth 
changes  without  imparting  a  pitch  angle  on  the  submarine. 
The  reader  should  note  that  M(t)  is  not  an  inertial 
reference  plane  variable,  and  represents  the  true  vertical 
rate  only  when  the  submarine  has  zero  pitch  and  roll 
angles.  Unfortunately,  z(t)  is  not  available  as  a  state 
variable  and  this  discrepancy  must  be  tolerated. 


With  the  output  variables  determined,  and  the  A  and  B 
matrices  calculated,  the  state  space  description  of  the 
submarine  model  is  oom  complete  and  takes  the  form 

x  <t>  »  A  x^t)  +  B  u(t> 

*<t>  *  C  x Ct >  ,  (2. 1 

where  y<t>  -  C  w<t>  r(t)  6<t)  3*. 


2.6  Pgrforiijinct  Requirements 


Two  performance  raquircmints  ar«  imposed  on  the  control lor 
design.  First,  ths  maximum  crossover  frequency  is  limited 


to  that  of  the  submar ii 


ability  to  respond  and,  second. 


the  steady  state  error  to  step  commands  and  step 
disturbances  is  to  be  zero. 


The  maximum  crossover  frequency  is  based  on  the  perturbed 
response  of  the  non-linear  model.  Inspection  of  figure  2.5 
reveals  an  average  settling  time  to  within  90%  of  steady 
state  of  60-65  seconds,  leading  to  a  desired  crossover 
frequency  of  0.10  radians/sec. 


Inspection  of  the  perturbed  response  of  all  four  models 
confirms  an  intuitive  feeling  that  the  settling  time 
increases  as  the  submarine  slows,  requiring  that  a  different 
performance  requirement  be  levied  on  each  design. 


The  value  of  0.10  rad/sec  was  declared  to  be  the  maximum 
for  the  fastest  model  (25  knots) ,  and  a  least  squares  curve 
was  fitted  to  the  approximate  time  constants  indicated  by 
the  plots.  A  linear  fit  proved  satisfactory,  giving  a 
correlation  coefficient  of  0.9983  and,  when  transformed  to 
yield  maximum  crossover  frequencies,  produced  the  results  in 


table  2.4 


Table  2.4 

Crossover  Frequencies  for  the  Linear  Models 

Model  Crossover  Frequency  (rad/sec) 

5  Knot  0.05 

10  Knot  0.06 

20  Knot  0.08 

25  Knot  0. 10 

The  zero  steady  state  error  requirement  is  met  by  placing 
integrators  in  each  of  the  three  input  channels.  The 
integrators  are  placed  at  the  input  to  the  plant  because 
this  is  where  the  error  signal  appears.  The  integrators 
will  then  become  a  part  of  the  compensator  which  is  before 
the  plant  in  the  feedback  loop.  As  will  be  seen,  the  use  of 
integral  control  does  not  prevent  the  maximum  crossover 
frequency  specification  from  being  met. 


2.7  Summary 


This  chapter  has  described  the  origins  of  the  submarine 
model,  briefly  described  the  process  of  producing  a  linear 
model,  and  presented  the  reasoning  by  which  the  output 
variables  were  chosen.  Additionally,  performance 
constraints  were  imposed  on  the  controller  design,  thus 
setting  the  stage  for  analysis  of  the  open  loop  dynamics  of 
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the  plant. 


The  next  chapter  will  analyze  the  plant  using  the  method 
of  modal  decomposition,  and  present  data  showing  the 
eigenstructure  of  the  submarine  linear  models. 


Chapter  3 


ANALYSIS  OF  THE  LINEAR  HODEL 


3.1  Introduction 


In  this  chapter,  the  eigenstructure  of  the  unaugmented  20 
knot  model  will  be  investigated.  As  previously  stated,  the 
20  knot  model  was  the  -first  to  be  developed  and  the  only  one 
to  have  a  detailed  analysis  performed  on  it. 

For  each  of  the  other  three  linear  models  (at  5,  10,  and 
25  knots)  analysis  was  limited  to  inspection  of  the  poles 
and  zeros,  and  producing  singular  value  plots  of  the  open 
loop  plant.  Several  consistency  checks  were  also  performed 
in  an  effort  to  avoid  any  numerical  errors. 

In  the  previous  chapter,  a  state  space  description  of  the 
submarine  linear  model  was  developed  in  the  form 

*(t>  «  A  x(t)  +  B  u ( t )  (3.1) 

£<t>  «  C  u  ( t ) .  (3.2) 

Prior  to  the  actual  analysis  of  the  linear  model,  it  will  be 
shown  that  the  order  of  the  model  can  be  reduced,  based  on 


the  zero  entries  present  in  the  A  and  B  matrices. 


The  eigenstructure  of  the  20  knot  model  is  revealed  by  the 
method  of  modal  analysis  C133.  This  method  involves  the 
implicit  application  of  similarity  transf ormations  to  the  A 
matrix,  reducing  it  to  a  diagonal  matrix  whose  entries  are 
the  poles  o-f  the  open  loop  system.  The  similarity 
transformations  provide  a  means  of  describing  the  state 
space  of  the  linear  system  in  separate  decoupled  modes, 
thereby  yielding  information  as  to  the  control  1 abi 1 i ty  and 
observability  of  the  system. 


This  information,  coupled  with  the  pole-zero  structure, 
provides  the  basis  and  validity  for  the  LQG/LTR  designs  of 
the  next  chapter. 


3.2  Reduction  of  the  Model 


Inspection  of  the  A  matrix  for  the  20  knot  model  (refer  to 
Appendix  A,  model  S20R5)  reveals  zero  entries  in  all 
elements  of  the  last  two  columns.  This  indicates  that  the 


| 


present  values  of  the  states  V  and  z  can  have  no  influence 
on  any  other  state.  In  other  words,  the  present  depth  and 
heading  angle  are  irrelevant  to  the  dynamic  response  of  the 
linear  system. 


m 
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Zeros  in  the  last  two  rows  uf  B  (Appendix  A)  further 
indicate  that  the  controls  exert  no  direct  influence  on  the 
derivatives  of  V  and  z. 

Therefore,  since  the  purpose  of  the  controller  is  not 
concerned  with  controlling  either  of  these  variables,  they 
are  removed  from  the  linear  model.  This  is  accomplished  by 
deleting  the  last  two  rows  and  columns  of  A,  and  the  last 
two  rows  of  B,  leaving  an  8  state  system. 

3.3  The  Natural  Modes  of  the  System 

Determination  of  the  natural  modes  of  the  20  knot  linear 
model  is  accomplished  by  diagonalizing  the  state  space 
description.  For  any  linear  dynamic  system  without  input  to 
output  direct  coupling,  we  have 

x(t>  =  A  x(t>  +  B  u  ( t )  (3.3 

*(t)  =  C  x(t).  (3.4 

Now  define  a  new  state  vector  £(t)  by  the  relation 

x  (t)  -  T  z  (t>  ,  (3.5 

where  T  is  an  as  yet  unspecified,  n  x  n,  invertible  matrix. 
Then  we  have 

T  z(t>  -  A  T  z(t>  +  B  u(t)  (3.6 

£(t)  «  C  T  z  (t) .  (3.7 

Multiplying  (3.6)  on  the  left  by  T  \  we  are  left  with 


z(t>  *  I  A  T  z(t)  +  I  B  u(t)  (3.8) 

y(t>  *  C  T  z  (t )  .  (3.9) 

I-f  T  is  such  that  T  *  A  T  is  diagonal,  than  the  state  vector 
z(t)  defines  a  new  state  space  in  which  the  modes  are 
decoupled  and  the  entries  of  T  1  A  T  are  the  eigenvalues  of 
A.  Let 


A  «  I  1  A 

T. 

(3.  10) 

Then 

we  have 

A  I  -  I  A. 

(3.11) 

Now, 

represent 

the  T  matrix  by 

its 

n  column  vectors  v^ ,  v^. 

■  ■  •  ? 

v^,  so  that  equation  3.11 

can 

be  expressed  as 

A  v.  =  v.x.  for  i  =  1  to  n  (A  =  U.)).  (3.12) 
-  -l  -l  i  i 

Thus  the  columns  of  T  are  the  eigenvectors  of  A  and  the 
diagonal izat ion  is  possible  only  when  the  A  matrix  has 
distinct  eigenvalues. 

The  T  matrix  is  called  the  modal  matrix .  with  each  column 
of  T  describing  the  motion  of  the  submarine  along  the 
coordinate  axes  of  the  state  vector  components  u,  v,  w,  p, 
q,  r,  ♦  ,  and  8,  for  a  particular  mode.  Since  every  possible 
dynamic  response  of  the  submarine  must  consist  of  a  linear 
combination  of  the  decoupled  modes,  analysis  of  the  columns 
of  T  should  provide  the  designer  with  very  useful 
information  regarding  the  submarine's  dynamic  response. 

Unfortunately,  it  is  not  a  simple  matter  to  visualize  or 


characterize  motions  in  8-space,  especially  when 


the 


magnitude  of  the  motions  are  complex  numbers  expressed  in 
various  units. 

Analysis  of  the  modes  of  the  system  is  accomplished  in  two 
steps.  First,  the  A  and  B  matrices  of  the  linear  model  are 
transformed  so  that  all  angular  components  of  the  state 
vector  will  have  units  of  degrees  or  degrees/sec,  as 
appropriate.  Details  of  the  transformations  applied  are  in 
appendix  Bl. 

Secondly,  the  columns  of  the  resulting  modal  matrix  are 
graphed  in  bar  chart  form  by  taking  absolute  values  of  each 
element  of  the  normalized  column  vectors.  The  loss  of  phase 
information  resulting  from  taking  absolute  values  is 
considered  to  be  of  less  importance  than  the  magnitude  of 
the  motion  along  a  particular  state  component. 

The  bar  charts  (figure  3.1)  have  a  vertical  scale  of  0  to 
100%,  reflecting  the  intent  of  displaying  the  relative 
magnitude  of  the  response,  and  the  eigenvalue  corresponding 
to  a  particular  mode  is  noted  below  each  graph.  Note  that 
the  two  columns  corresponding  to  a  complex  eigenvalue  do  not 
produce  unique  modes  when  graphed  in  this  manner.  Complex, 
or  oscillatory  modes  are  indicated  by  the  labeling  and  the 
presence  of  a  complex  eigenvalue. 

It  should  be  noted  at  this  point  that  the  modes  shown 


represent  the  natural  modes  of  a  linear  model  of  a  submarine 


proceeding  downward,  in  a  port  turn,  at  20  knots.  (Recall 
that  the  nominal  point  specifies  5  degrees  of  rudder  and  1 
degree  on  the  planes).  Thus,  although  the  bar  charts 
provide  a  convenient  means  of  displaying  the  modes  of  a 
linear  system,  in  this  particular  case  they  are  displaying 
the  modes  of  a  system  linearized  about  a  nominal  point  that 
serves  to  further  obscure  any  physical  interpretation.  For 
this  reason,  interpretation  of  the  modes  is  limited  to  the 
following  observations: 

1.  All  eigenvalues,  and  hence  the  open  loop  poles,  are  in 
the  left  half  plane. 

2.  Mode  1  appears  to  represent  the  steady  state  response 
of  the  submarine  in  a  diving  turn  with  constant 
control  surface  deflections. 

3.  Modes  1,  2,  and  3  are  relatively  lightly  damped 

compared  to  the  other  modes.  They  also  have  less  roll 
response  and  higher  pitch  response  than  the  others, 
and  tend  to  indicate  that  the  sail  of  the  submarine 
acts  as  a  damper  in  side  to  side  motions  of  the 
submarine. 

4.  Modes  5  and  6  represent  a  nearly  pure  roll  mode  in 
that  almost  the  entire  response  is  concentrated  in 
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roll  and  roll  rat*.  The  inability  to  directly  control 


roll  and  roll  rate  is  reflected  in  the  oscillatory 
behavior  of  these  modes. 

5.  Comparison  of  the  magnitudes  for  roll,  pitch,  and  yaw 
rates  for  modes  4  and  8  appears  to  indicate  that  as 
more  yaw  rate  is  introduced,  some  of  the  roll  response 
is  converted  to  pitch  response.  This  is  indicative  of 
the  cross  coupling  that  exists  between  the  stern 
planes  and  the  rudder  when  the  submarine  is  rolling 
and  turning  at  the  same  time. 

The  eigenvalues  and  modal  matrices  for  all  models  are 
contained  in  appendices  B2  through  B5. 


3.4  Controllability  and  Observability 

The  LQG/LTR  design  methodology  provides,  under  certain 
conditions,  guarantees  to  the  control  system  designer 
concerning  gain  and  phase  margins  C8J.  It  should  be  noted 
that  the  guarantees  are  not  absolute,  and  exist  only  in  the 
sense  that  the  LTR  phase  of  the  design  process  recovers  a 
loop  shape  that  approaches  that  of  the  Kalman  Filter  design, 
which  does  have  the  robustness  guarantees. 


The  linear  system  must  meet  certain  requirements,  however. 


before  the  validity  of  the  results  are  assured. 
Specifically,  the  system  must  have  no  unstable  modes  which 
are  not  controllable,  and  no  unstable  modes  which  are  not 
observable.  If  controllability  and  observability  of  the 
system  can  be  established,  then  the  weaker  conditions  of 
stabi 1 izabi 1 i ty  and  detectability  are  assured. 

When  a  linear  system  can  be  diagonalized,  as  in  the 
present  case,  determination  of  control labi 1 i ty  and 
observability  becomes  a  very  simple  matter.  Additionally, 
even  when  a  system  does  not  meet  the  criteria  for  being 
controllable  or  observable,  the  weaker  conditions  of 
stabi 1 izabi 1 ity  and  detectability  are  also  readily 
determined. 

Recall  that  the  new  state  vector  z_  (equation  3.5)  defines 
a  non-physical  state  space  in  which  the  natural  modes  of  the 
system  are  decoupled.  Thus  a  particular  row  of  the  T  1  B 
matrix  links  the  input  vector  u  to  a  particular  mode.  Each 
element  in  the  row  then  links  a  particular  input  (in  our 
case  a  control  surface)  to  a  mode.  Therefore,  a  zero  entry 
in  the  (i,j)  position  of  the  T  1  B  matrix  would  indicate 
that  the  ith  mode  is  not  controllable  by  the  jth  input. 

Similarly,  the  matrix  C  T  in  equation  3.9  indicates 
whether  a  particular  mode  is  visible  in  the  output. 


a 


i 

& 


The  20  knot  linear  model  has  non-zero  entries  in  the  T  B 
and  C  T  matrices  (refer  to  appendix  B4>  and  is  therefore 
both  controllable  and  observable,  providing  us  with 
assurances  that  the  LQG  design  method  to  be  applied  in  the 
next  chapter  will  have  no  difficulties. 


3.5  Poles,  Zeros,  and  Singular  Values 


The  poles  of  the  open  loop  plant  are  the  eigenvalues  of 
the  A  matrix  presented  in  the  previous  section.  As  was 
indicated,  the  plant  is  open  loop  stable,  having  all  left 
half  plane  poles. 


Multivariable  transmission  zeros  appearing  in  the  right 
half  plane  may  present  difficulties  in  an  LQG/LTR  design 
C81.  Fortunately,  the  20  knot  linear  model  has  no 
non-minimum  phase  zeros,  nor  do  any  of  the  other  models 
(refer  to  Appendix  B6) . 


The  open  loop  poles  for  all  four  models  are  presented  with 
the  modal  information  in  appendices  B2  through  B5,  and  the 
transmission  zeros  for  each  model  in  Appendix  B6. 


Singular  value  plots  of  the  open  loop  models  are  shown  in 
figure  3.2.  It  will  be  noticed  that  integral  control  is  yet 
to  be  added,  and  that  satisfactory  command  following 
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performance  will  indeed  require  the  addition  of  integrators 
in  the  loop  transfer  -function  matrix.  The  increased  e-f-fect 
o-f  the  control  surfaces  as  speed  increases  can  be  seen  by 
noting  the  increased  dc  gain  for  each  sucessive  model. 

As  a  further  indication  of  the  effect  of  the  control 
surfaces  on  the  outputs,  the  dc  gains  of  the  open  loop 
transfer  function  matrix  for  the  20  knot  model  are  listed  in 
table  3.1.  By  reading  across  for  each  output  variable,  the 
relative  effect  of  the  control  surfaces  can  be  determined. 
The  results  are  as  expected,  with  the  sail  planes  having  the 
most  effect  on  w,  the  rudder  on  r ,  and  the  stern  planes  on 
©.  The  gains  for  w  are  significantly  higher  than  for  the 
other  two  outputs,  indicating  a  need  for  scaling,  which  is 
done  in  the  next  chapter. 

Table  3.1.  Input  to  Output  Coupling 


ds 

db 

dr 

w 

11  db 

19  db 

16  db 

r 

-56  db 

-49  db 

-18  db 

e 

19  db 

-8.1  db 

-3.0  db 
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This  chapter  on  analysis  of  the  linear  model  would  be 
incomplete  without  demonstrating  the  extreme  usefulness  of 
modal  analysis  in  detecting  inconsistencies  in  a  model. 

The  initial  linear  model  of  the  submarine  on  the  CSDL 
computer  contained  an  11th  state  known  as  n.  The  purpose  of 
this  state  was  to  include  the  propulsion  plant  dynamics  and 
was  defined  as  u (actual ) /u (commanded) . 

Modal  analysis  of  the  model  with  *  included  resulted  in 
eight  modes  essentially  the  same  as  those  depicted  in  figure 
3. 1 ,  and  a  ninth  mode  with  all  response  concentrated  in  u 
and  n-  Furthermore,  this  ninth  mode  had  a  near-zero 
eigenvalue  and  zeros  in  the  corresponding  row  of  the  T  1  B 
matrix . 

From  the  previous  discussion  of  modal  analysis,  we  can  see 
that  this  situation  is  indicative  of  an  uncontrollable  and 
unstable  mode,  thus  rendering  the  LQ6  design  methodology 
invalid.  If,  however,  the  near — zero  eigenvalue  is  taken  to 
be  the  digital  representation  of  a  true  zero  eigenvalue,  and 
the  reasoning  is  applied  that  the  actual  submarine  must  be 
stable  in  forward  velocity,  then  the  facts  imply  an 
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inconsistency  in  the  model. 

To  further  refine  this  conclusion,  both  the  linear  and  the 
nonlinear  equations  of  motion  were  integrated  over  time, 
using  the  scaled  magnitudes  of  the  eigenvector  corresponding 


to  the  unstable 

mode.  The  nonlinear  plot 

decayed 

back 

to 

the  local 

equi librium 

point  whereas 

the  linear 

plot 

di spl ayed 

an 

unf orced 

response  that 

remained 

at 

the 

perturbed  values. 

This  last  information  clearly  indicated  a  dependency  among 
the  states  of  the  linear  model,  and  the  problem  was  then 
quickly  traced  to  the  treatment  of  the  propulsion  plant  as  a 
constant  rpm  source  for  the  propeller.  Thus  the  commanded  u 
was  being  taken  as  constant  for  any  run,  and  n  represented  a 
constant  multiple  of  u,  producing  the  dependent  state. 

Application  of  the  chain  rule  to  the  system  of  11 
equations  in  11  unknowns  (used  to  calculate  the  A  and  B 
matrices)  reduced  the  order  of  the  system  to  10,  and 
eliminated  the  unstable  state. 

It  is  hoped  that  this  example  may  serve  to  aid  the  reader 
in  the  methodical  elimination  of  similar  problems  in  a 
linear  model. 
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3.7  Summary 

This  chapter  has  described  the  technique  of  modal  analysis 
and  its  use  in  determining  the  eigenstructure  and  modal 
composition  of  the  state  space  description  of  a  linear 
model . 

Application  of  this  technique  has  enabled  us  to  establish 
the  prerequisi tes  required  to  pursue  the  LQG/LTR  design 
procedure  to  be  covered  in  the  next  chapter.  These  are  that 
the  open  loop  linear  model  be  detectable  and  stabi 1 izable, 
and  that  the  location  of  any  non-minimum  phase  zeros  be 
determined. 

Finally,  an  example  of  the  usefulness  of  modal  analysis  in 
tracking  down  an  error  in  the  linearization  of  the  nonlinear 


model  was  presented 


Chapter  4 

LINEAR  MIHO  0ESIGN8 

4.1  Introduction 

This  chapter  presents  the  linear  portion  of  the  design 
process,  utilizing  the  LQG/LTR  design  methodology  as  the 
NIMO  design  tool. 

The  design  section  begins  with  an  overview  of  the  complete 
step  by  step  LQG/LTR  procedure,  as  tailored  to  conform  to 
the  main  focus  of  this  thesis. 

Augmentation  of  the  submarine  model  dynamics  with 
integrators  will  then  be  discussed,  along  with  a  comparison 
of  open  loop  singular  value  plots  for  the  augmented  and 
unaugmented  20  knot  model . 

Finally,  the  methods  used  to  ensure  the  identical  behavior 
of  the  singular  values  of  the  open  loop  plant  at  both  the 
high  and  low  ends  of  the  frequency  spectrum  will  be 


Testing  of  the  compensator  begins  with  a  comparison  of  the 
step  responses  of  both  the  linear  and  nonlinear  models,  and 
concludes  with  simulations  of  evasive  maneuvers  and  control 
surface  saturation  tests. 

All  parameters  used  in  the  design  process  and  many  of  the 
design  products  for  all  four  models  are  contained  in  the 
appendi ces. 


4.2  The  LQG/LTR  Design  Methodology 

The  multivariable  LQG/LTR  design  methodology  consists  of  a 
four  step  process  £83. 

Step  one  involves  the  development  of  a  low  frequency  model 
of  the  nominal  plant  and  determination  of  the  uncertainties 
present  in  the  model.  It  is  assumed  that  the  frequency 
range  of  interest  for  purposes  of  command  following  and 
disturbance  rejection  is  at  low  frequencies. 

The  uncertainty  in  the  nominal  model,  including  sensor 
noise,  unmodeled  dynamics,  and  certain  actuator  dynamics,  is 
similarly  assumed  to  be  concentrated  at  high  frequencies. 
The  determination  of  what  constitutes  high  and  low 
frequencies  consists  of  fixing  the  maximum  allowable 


crossover  frequency  for  the  nominal  design 


In  this  thesis,  the  actual  linear  time  invariant  plant  and 


the  low  frequency  nominal  model  are  taken  to  be  identical, 
and  no  attempt  is  made  to  establish  the  modeling 
uncertainty.  Therefore,  step  one  is  restricted  to  producing 
the  linear  model,  as  was  discussed  in  chapter  2,  and 
determination  of  the  maximum  allowable  crossover  frequency 
for  each  model,  also  discussed  in  chapter  2. 


Step  two  of  the  design  process  establishes  the  low 
frequency  performance  requirements.  The  frequency  domain 
block  representation  of  the  compensated  plant  is  shown  in 
figure  4.1. 


Figure  4.1.  Block  Diagram  of  a  MIMO  Compensated  Plant 


r  <s> 
e  <s> 
u  (s) 

(s) 

P 

d  (s) 
K(s) 
6  (s) 


reference  signal  or  command  input  vector 
error  signal  vector 
control  input  vector  to  the  plant 
output  vector  of  the  plant 

disturbance  vector  reflected  at  the  plant  output 
compensator  transfer  function  matrix 
(augmented)  plant  transfer  matrix 
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The  transfer  matrix  G(s)  is  assumed  to  contain  the  nominal 

low  frequency  model  of  the  submarine  §p<*>  plus  any 

augmented  dynamics  Ga<s)  and  is  termed  the  nominal  design 

model.  Thus  G(s)  *  G  (s)  G  (s>.  To  determine  the 
—  -p  — a 

requirements  to  be  imposed  on  the  compensator  transfer 
matrix  K(s> ,  the  overall  transfer  function  of  the  closed 
loop  system  is  calculated 

^(s)-Ci+GCs)K(s) 3-1d (s) +  C  I^-G  <  s ) K  <  s  > 3~*G(s)K(s)r(s).  (4.1) 

Good  command  following  requires  that  £(s)  ~  r(s>  and  good 
disturbance  rejection  will  result  if  the  contribution  of 
d(s)  can  be  kept  small.  Here,  s  -  jw  is  in  the  low  end  of 
the  frequency  spectrum. 

Both  of  these  requirements  can  be  met  if  the  minimum 
singular  values  of  the  matrix  product  G(s)K(s)  are  large 
with  respect  to  unity  at  frequencies  below  crossover. 

Similarly,  the  response  of  the  outputs  to  high  frequency 
effects  can  be  minimized  if  the  maximum  singular  values  of 
G(s)K(s)  are  small  compared  with  unity  at  frequencies  above 
crossover . 

The  combined  effect  of  steps  1  and  2  is  to  place  barriers 
on  the  singular  value  plot  for  6(jw)K(jw)  as  indicated  in 
figure  4.2.  The  high  frequency  barrier  places  a  robustness 
requirement  on  the  compensator  and  the  low  frequency  barrier 


is  a  result  of  the  command  following  and  disturbance 


rejection  requirements. 


20  log  <7 


The  remainder  of  the  design  process  is  concerned  with 
determining  the  transfer  matrix  K(s)  so  that  the  singular 
values  of  the  open  loop  transfer  matrix  G(jw)K(jw)  will  have 
the  shapes  shown  in  figure  4.2.  This  process  is  known  as 
loop  shaping  and  is  at  the  heart  of  the  LQG/LTR  design 
procedure. 


Initially,  the  Kalman  Filter  methodology  is  applied  to  the 
state  space  description  of  the  nominal  design  model  (nominal 
model  plus  augmentation)  to  produce  a  transfer  matrix  G^-Cs) 
that  has  the  desired  loop  shapes.  The  Kalman  Filter  theory 
is  applied  in  a  very  specific  manner  and  should  not  be 
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con-fused  with  optimal  state  estimation  applications 


To  solve  the  Kalman  Filter  problem,  the  nominal  design 
model  state  space  description, 

x(t)  -  A  x(t)  B  u(t)  (4.2) 

%<t)  »  C  x <t) ,  (4.3) 

is  re-formulated  to  produce  the  somewhat  modified  dynamics 

x(t)  «  A  x(t>  *  L  £(t)  (4.4) 

£(t)  *  C  x(t)  +  0(t).  (4.5) 

where  Y(t)  represents  the  process  white  noise, 
and  0(t)  represents  the  measurement  white  noise. 

The  covariance  of  these  two  Gaussian  signals  is 

covCX  (t)  li(T)  3  *  U(t-t>,  (4.6) 

and  covte(t)  |€>(t)  3  -  *I&(t-T>  (/i  >  0)  .  (4.7) 


The  design  parameters  a  and  L  are  then  used  to  produce  the 
desired  loop  shapes  of  the  transfer  matrix  G^p(s)  where 

G^s)  »  C(sI^-A)~1H,  and  (4.8) 

H  «  (1/V>>  E  C'.  (4.9) 

The  matrix  E  is  the  solution  to  the  Filter  Algebraic 
Riccati  Equation  (FARE) 

0-AE+EA'  ♦  L  L'  -  (1/V>>  E  C'C  E.  (4. IB) 

Fortunately,  the  transfer  matrix  G^p(s)  can  be  approximated 
quite  readily,  providing  that  the  parameter  ji  <<  1,  and  the 

FARE  need  not  be  evaluated  for  each  choice  of  jl  and  L.  It 
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can  bs  shown  ZB1  that 

q^CG^pUw)]  *  <1/V>>  q^  CGpgL  ( j»)  3 ,  -for  ft  <<  1,  (4.11) 

where  §p0L<j<*>  ■  C (sI^-A) _1L.  (4.12) 

Thus  tha  L  matrix  can  ba  chosen  in  such  a  way  as  to 
produce  tha  desired  loop  shape  and  /l  used  to  adjust  tha 
singular  values  up  and  down  to  meet  the  required  crossover 
■frequency  criteria. 


Providing  that  CA,LJ  is 
detectable,  any  choice  of 
guaranteed  properties  for 


stabi 1 i zabl 
ft.  and  L  will 


e,  and 
produce 


that  CA,C3  is 
the  following 


1.  Closed  loop  stable 

2.  Robustness i 

<yicI+®KF( J**>  3  3k  i  (4.13) 

ajCI+G^-p-1  (jw)  3  *  1/2  (4.14) 


3.  Infinite  upward  gain  margin 

4.  6  db  downward  gain  margin 

5.  Phase  margins  of  ±  60* 


The  final  step  in  the  LQ6/LTR  design  procedure  involves 
the  "recovery"  of  the  loop  shapes  of  G^p(s)  by  the 
compensated  plant  transfer  matrix  G(s)K(s>.  This  is 
accomplished  by  solving  the  Control  Algebraic  Riccati 
Equation  (CARE) 
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0  »  -  K  A  -  A'K  -  qC'C  K  B  B'K  j  q  >  0  (4.15) 
using  ths  dssign  parameter  q,  and  defining  the  control  gain 
•atrix  G  by 

G  *  B'K.  (4.16) 

For  the  solution  to  the  CARE  to  be  valid,  mb  must  have 
CA,B1  stabilizable  and  CA,C3  detectable.  Furthermore,  the 
nominal  design  plant  must  not  have  any  non-minimum  phase 
zeros. 


The  control  gain  matrix  G  and  the  filter  gain  matrix  H 
(4.9),  when  calculated  using  the  above  procedures,  define  a 
special  type  of  compensator  known  as  a  Model  Based 
Compensator  (figure  4.3).  This  type  of  compensator  differs 
from  other  LQG/LTR  compensators  only  in  the  manner  in  which 
the  control  and  filter  gain  matrices  are  calculated. 


The  state  space  description  o-f  the  (model  based) 


compensator  is  then  given  by 


z  ( t )  -  <  A-B  G-H  C)  z ( t )  -  H  e<t) 


u(t)  *  -G  z  <t> , 


<4. 17> 


(4.  18) 


It  can  be  shown  that  the  singular  values  of  G(5)KM__(s)  at 

—  -rloU 

frequencies  below  crossover  converge  to  those  of  G^pls)  as 
the  design  parameter  q  -»  co,  providing  that  the  plant  is 
minimum  phase  C81.  At  frequencies  above  crossover,  an  extra 
pole  of  roll-off  is  produced  by  the  recovery  phase,  further 
enhancing  the  high  frequency  robustness.  Thus  the  loop 
shape  of  G^ts)  is  recovered  and  the  resulting  controller 
will  have  the  desired  performance  characteristics  described 
ear 1 i er . 


ft 


i 


v| 
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4.3  Augmented  Dynamics 


The  dynamics  of  the  submarine  model  are  augmented  by 
placing  integrators  in  each  of  the  three  control  channels. 
This  is  accomplished  by  defining  an  augmentation  plant  G^  (s) 
whose  state  space  description  is  simply 


A  »  0,  B  -  C  -  I, 

-a  —  —a  —a  — 


with  all  matrices  being  3x3. 
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Figure  4.4.  Inclusion  of  the  Augmented  Dynamics 


Then  the  nominal  low  frequency  model  and  the  augmentation 

dynamics  are  combined  by  per-forming  a  state  space 

multiplication,  forming  G(s>  as  shown  in  figure  4.4.  Note 

that  the  physical  input  to  the  plant  is  now  labeled  u  (s)  to 

~P 

distinguish  it  from  the  output  of  the  compensator  u(s) . 
Although  the  augmentation  dynamics  G^(s)  will  eventually  be 
lumped  with  the  compensator,  they  are  kept  separate  until 
the  LQG/LTR  design  procedure  is  complete.  Figure  4.5  shows 
a  comparison  of  the  singular  values  of  the  unaugmented  and 
augmented  20  knot  model.  As  can  be  seen,  adding  the 
integrators  at  the  plant  input  produces  a  60  db  gain 
increase  at  .001  rad/sec. 


2.4  Kalman  Filter  Design 

In  the  overview  of  the  LQG/LTR  design  methodology,  it  was 
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stated  that  the  singular  values  of  the  Kalman  Filter 
transfer  function  matrix  are  closely  approximated  by  the 
singular  values  of  ( 1/V>.>  G^g^  ( J»> .  Furthermore,  Gp0L<jw> 
can  be  easily  calculated  for  each  choice  of  L,  simplifying 
the  iterative  nature  of  the  design  process. 


To  satisfy  the  loop  shaping  requirements  represented  by 
the  low  and  high  frequency  barriers  shown  in  figure  4.2,  it 
is  helpful  if  the  maximum  and  minimum  singular  values  of 


are  identical  at  low  and  high  frequencies. 
Therefore,  the  choice  of  L  is  based  on  this  philosophy. 


Recall  that  G(s)  =  G  <s)G  <s>  (section  4.3>  and  define  a 

-a  -p 

state  space  representation  for  G(s)  by  the  relation 


G(s)  *  C ( s T-A )  B. 

Then  the  matrices  A  and  C  can  be  partitioned  into 


C  -  C  0  C  J, 

—  —  -p 


and  similarly  for  s£-A  and  its  inverse 


(4. 19) 


(4.20) 


sI-A 


si  0 


-B_  si -A 


(si -A) 


(si -A  > _1B  /s  (sI-A)  1 


(4.21) 


At  low  frequencies,  si -A  »  -A  and  (sI-A_ >~*  *  -A 

- p  “P  - P  -p 

Recall  that  A  *  exists  since  A  has  distinct  and  non-zero 
“P  “P 

eigenvalues.  Me  now  partition  the  L  matrix  into  Lj  and 
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-v1 


(4.22) 


and  form  G__.  (s)  for  low  frequencies 
— FOL  _ 


§__,<*)  *  C(sl^-A)  _1L  »  C0  Cl 

-A  -1B  /i 
_~P  ~P 

*  -C  A  _1B  L. /*  -  C  A  “1L_. 

-P  “P  -P-1  -P-P  “2 


From  this,  it  is  saan  that  wa  can  match  tha  singular 
valuas  at  low  (re quancy  by  choosing 


L,  -  -<C  A  *B  )  1, 

-1  -p-p  ~p 


(4.23) 


i-f  it  exists,  and  leaving  unspecified.  Then,  as  u  -»  0, 


GrQL  ( J«»)  ~  I/j«  ♦  M, 


(4.24) 


where  M  is  a  3x3  constant  matrix.  Although  M  will  have  an 
impact  at  -frequencies  approaching  crossover,  the  first  term 
will  dominate  for  <■>  small  enough. 

At  high  frequencies,  sI-A  ~  sl^  and  (s£-A  )  *  ~  I/s.  Thus, 


at  high  frequencies,  we  have 


W> 


1 

pJ  B  /ss 
“P 


i/s  l2 


(4.25) 


C  B  L,/s* 

-p-p-1 


W 


Since  1/s  >  l/s*  as  s  -»  os,  than  tha  second  term  will 
dominate  the  maximum  singular  value  at  high  frequencies. 
Thus  wa  can  match  tha  singular  valuas  at  high  frequency  by 
choosing 
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-2  -p  -p-p 


(4.26) 


so  that,  as  *>  -»  oo. 


G__.  < jm)  •  -C  B  L,/w*  +  I/-. 

-FOL  -p-p— 1  — 


(4.27) 


The  required  invtrsn  can  ba  shown  to  axist  for  tha  models 


usad  via  nuaarical  aaana.  Additionally,  whenever  C  is  such 

“P 

that  distinct  states  are  extracted  from  the  state  vector, 
the  matrix  C  C  ‘  ii  always  diagonal  and  hence  invertible. 


The  above  method  of  constructing  the  L  matrix  provides  the 
designer  with  a  guarantee  of  identical  behavior  of  the 
singular  values  of  the  Kalman  Filter  transfer  matrix  for  low 
and  high  frequencies.  Unfortunately,  the  shape  of  the 
singular  values  at  frequencies  near  crossover  is  not 
controlled,  and  significant  differences  may  exist  between 
the  maximum  and  minimum  singular  values  in  the  crossover 
region.  Figure  4.6  is  a  plot  of  the  singular  values  of 
GpQL(s)  for  the  L  matrix  constructed  as  in  equations  4.23 
and  4.26.  Although  the  singular  values  match  at  high  and  low 
frequencies,  a  large  separation  exists  between  the  maximum 
and  minimum  crossover  frequencies. 


To  control  the  separation  of  the  singular  values  at 
frequencies  near  crossover,  it  is  necessary  to  investigate 
the  magnitudes  of  the  elements  of  the  matrix  transfer 
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function  Gpg^(j<i»)  over  a  frequency  range  including  the 
crossover  frequency,  and  apply  scaling  to  the  matrix  to 
achieve  a  tight  crossover  pattern.  The  magnitudes  of  the 
elements  of  the  matrix  transfer  function  are  plotted  in 
figure  4.7  for  each  input  to  output  function.  Each  of  the 
plots  represents  the  coupling  from  the  three  controls  to  a 
single  output.  The  curves  are  fairly  smooth  except  for 
figure  4.7(a),  which  shows  the  coupling  from  the  controls 
ds,  db,  and  dr  to  the  output  w.  The  desired  crossover  is  at 
u  =  0.100,  or  log  fa>  =  -1.00,  and  figure  4.7(a)  displays  both 
an  increase  in  coupling  from  dr  and  a  decrease  in  coupling 
from  db  in  this  area.  Referring  back  to  table  3.1,  it  was 
observed  that  the  dc  gains  for  w  were  significantly  higher 
than  for  r  and  0.  Both  of  these  facts  tend  to  indicate  that 
the  Cp  matrix  entry  corresponding  to  w  should  be  scaled 
down.  A  comparison  of  the  magnitudes  of  the  transfer 
matrices  in  this  manner  for  all  four  submarine  models 
resulted  in  the  following  choice  for  the  C  matrix. 


C 

-P 


00.1  00000 
00000100 
00000001 


(4.28) 


The  vast  improvement  in  the  behavior  of  the  singular 
values  near  crossover  can  be  seen  in  figure  4.10(a),  a  plot 
of  the  singular  values  of  6^ 
using  the  new  C  matrix. 


(s>  for  the  20  knot  model 


(b)  Kalman  -filter  open  loop  G(kf) 


I  mm  »m»t« 

(c)  Recovered  open  loop  transfer  function 
Figure  4.9.  Design  Summary  for  the  10  Knot  Model 
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(c)  Recovered  open  loop  transfer  function 
Figure  4.11.  Design  Summary  for  the  25  Knot  Model 
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submarine  models  did  not  have  any  low  frequency  transmission 


zeros.  A  value  of  q  -  100.0  was  used  for  all  models, 
producing  the  control  gain  matrices  in  appendix  C2. 

The  complete  design  sequence  is  summarized  in  figures  4.8 
through  4.11,  the  singular  value  plots  of  6^.^(s),  6^„(s)  , 
and  G(s)K(s)  for  each  of  the  four  models.  The  maximum  and 
minimum  crossover  frequencies  for  the  final  loop  transfer 
function  matrix  C6(s)K(s)l  are  summarized  in  table  4.2. 

Table  4.2.  Final  Crossover  Frequencies  for  G(s)K(s) 


Model 

u 

mm 

w 

max 

S5R5 

.016 

.050 

S10R5 

.020 

.068 

S20R5 

.029 

.078 

S25RS 

.024 

.079 

4. A  The  Closed  Loop  System 

Prior  to  beginning  the  testing  phase  of  the  compensator 
design,  properties  of  the  closed  loop  plant  will  be 
investigated. 

The  state  space  representation  of  the  closed  loop  system 
can  be  written  down  by  inspecting  figures  4.1  and  4.3.  Let 


x^(t )  represent  the  state  of  the  nominal  design  model  and 


z(t)  the  state  of  the  compensator, 
system  can  be  described  by 


X  <t) 

A 

-B  G 

z(t) 

H  C 

A-B  G-H  C 

£<t> 


x  (t) 
z  (t) 


Then  the  closed  loop 


■* 

* 

x  (t) 

0 

+ 

z  <t> 

-H  C 

_  — 

. 

<4.29) 


The  poles  and  zeros  for  the  closed  loop  system  are 
contained  in  appendix  C3.  All  poles  are  in  the  left-half 
plane  so  that  the  system  is  in  fact  closed  loop  stable. 


By  setting  d(s)  equal  to  zero  in  equation  4.1,  it  is  seen 
that  the  singular  values  of  the  closed  loop  plant  should  be 
very  near  unity  (zero  db  on  a  log  plot)  from  dc  up  to  the 
crossover  frequency,  and  then  fall  off  above  crossover. 
This  is  depicted  in  figures  4.12  and  4.13,  the  closed  loop 
singular  value  plots  for  all  four  models. 


4.7  Testing  of  the  Compensated  Submarine  Model 

Testing  of  the  LQG/LTR  model  based  compensator  design  was 
accomplished  by  providing  the  computer  program  with  a  data 
file  containing  time  sequenced  command  inputs  and  then 
integrating  either  the  linear  or  nonlinear  equations  of 


log  omega 


(b)  Tha  IB  knot  modal 
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motion.  The  command  inputs  -for  each  test  are  superimposed 
on  the  output  variable  time  histories  to  indicate  the  time 
that  the  command  Mas  applied  and  also  to  provide  an  implicit 
error  history. 


Results  are  shown  for  the  20  knot  model  only  as  it  enjoyed 
the  most  exhaustive  testing.  Check  runs  were  made  for  the 
other  models  and  displayed  similar  results. 


4.8  ComDari son  of  Linear  and  Nonlinear  Simulations 


The  initial  testing  consisted  of  a  comparison  of  the 
responses  of  the  linear  and  nonlinear  equations  of  motion  to 
step  inputs.  Figure  4.14  is  representat i ve  of  this  stage  of 
testing  and  shows  the  time  histories  of  the  three  output 
variables  (w,  r,  and  0)  and  the  three  control  surfaces  when 
the  submarine  model  is  subjected  to  step  inputs  in  w  and  8. 
For  this  plot,  the  initial  values  of  the  command  inputs  were 
set  to  the  nominal  values  for  the  20  knot  model ,  with  the 
step  changes  being  applied  at  t  s  50  seconds,  and  removed  at 
t  *  100  seconds. 


The  reader  should  note  the  coupling  that  exists  between 
the  stern  and  fairwater  planes  as  they  deflect  in  opposite 
directions  (as  they  should)  in  response  to  the  pitch  command 


(figure  4.14  j  and  1).  Additionally,  the  pitch  command 
completely  dominates  the  response  initially  and  the  command 
in  w  has  no  effect  until  the  pitch  error  decreases.  This 
should  not  be  the  case  if  the  inertial  variable  z  was  used 
in  place  of  w,  although  any  demonstration  of  this 
supposition  will  have  to  wait  until  the  computer  model  is 
changed. 

Figure  4.14  also  displays  the  effects  of  the  cross 
coupling  that  exists  between  yaw  rate  and  pitch,  as 
indicated  by  the  variation  in  yaw  rate  (r)  (figure  4.14  (e) 
and  (f))  and  the  rudder  deflection  (figure  4.14  (m)  and  (n) ) 
at  the  time  when  the  step  command  is  applied  in  6  (figure 
4.14  (g)  and  (h) ) .  This  is  a  result  of  both  the 
unconventional  operating  point  that  the  model  was  linearized 
about  and  the  non-zero  roll  angle.  Recall  that  the  nominal 
point  reflects  a  diving  turn  so  that  a  roll  angle  is  present 
for  this  maneuver. 

The  main  objective  of  this  test  was  to  determine  the 
extent  of  agreement  between  the  linear  and  nonlinear 
compensated  models  in  order  to  establish  the  validity  of 
both  the  compensator  design  and  the  computer  software  being 
used  for  the  simulations.  The  extremely  close  agreement 
between  the  time  histories  for  the  linear  and  nonlinear 
models  displayed  by  figure  4.14  thus  establishes  the 


validity  and  permits  -further.  All  -further  simulations  are 
per-formed  using  the  nonlinear  computer  model. 

4.9  Further  Simulation  o-f  the  Nonlinear  Model 

Having  established  the  validity  o-f  the  compensator  design 
by  successful  comparison  o-f  the  linear  and  nonlinear 
simulations,  the  next  step  is  to  examine  such  -factors  as 
symmetry,  multiple  commands,  and  control  surface 
saturati on. 

The  degree  of  symmetry  present  in  the  submarine  model  was 
evaluated  by  commanding  a  turn  first  in  one  direction  and 
then  in  the  opposite  direction.  The  results  are  shown  in 
figure  4.15  and  reveal  that  even  with  the  unconventional 
operating  point  used  in  producing  the  linear  model, 
remarkable  symmetry  exists  in  the  expected  state  variables. 
The  roll  angle  ($)  time  history  is  included  in  this  figure 
to  display  the  submarine's  tendency  to  lean  into  a  turn.  A 
rather  curious  phenomenon  is  also  displayed  in  that  the 
submarine  pitches  up  initially  when  commanded  to  turn.  This 
trait  reportedly  exists  in  an  actual  submarine  and  was  also 
observed  on  the  real  time  simulation  facility  mentioned  in 
the  introduction  to  chapter  2.  Finally,  notice  that  the 
performance  of  the  controller  begins  to  degrade  at  about  135 
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(b)  Forward  velocity  u(t>  -for  the  nonlinear  model 


(c)  Vertical  rate  w(t)  -for  the  linear  model 


0  20  U0  60  80  100  120  1U0  160  180 

TIME 

(d)  Vertical  rate  w(t)  for  the  nonlinear  model 
Figure  4.14.  Comparison  of  the  Linear  and  Nonlinear  Simulations 
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seconds  as  the  forward  velocity  reaches  33 7.  below  the 
nominal  point,  indicating  that  gain  scheduling  is 
necessary. 

The  multiple  command  test  was  accomplished  by  simulating 
an  evasive  maneuver.  The  command  inputs  are  applied  at 
different  times  and  simulate  first  a  diving  turn  followed  by 
a  recovery  of  the  initial  conditions.  The  results  of  this 
test,  shown  in  figure  4.16,  display  the  ability  of  the 
controller  to  respond  to  commands  in  both  yaw  rate  and 
pitch.  Again  we  see  that  the  pitch  command  has  significant 
effect  on  the  yaw  rate  due  to  the  cross  coupling.  Also 
displayed  is  the  inability  of  the  sail  planes  to  control  w 
in  the  presence  of  pitch  commands.  The  sail  planes  do, 
however,  aid  the  stern  planes  in  the  initial  response  to  a 
pitch  command  by  deflecting  in  the  opposite  direction,  as  is 
desired.  Note  the  loss  of  control  in  w  at  about  75  seconds 
when  the  sail  planes  saturate. 

To  fully  display  the  behavior  of  the  compensated  nonlinear 
model  during  periods  of  control  surface  saturation,  the 
reference  inputs  for  w  and  9  used  in  the  previous  test  were 
increased  significantly  in  magnitude  and  the  simulation 
repeated.  The  results  are  presented  in  figure  4.17  and 
display  a  loss  of  control  during  saturation  along  with 
periods  of  lightly  damped  oscillatory  response.  Although 
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Figure  4. 15.  Yaw  Rate  Symmetry  Test 
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(a)  Forward  velocity  u(t) 
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(b)  Vertical  rate  w(t) 


the  command  inputs  are  removed  at  150  seconds*  the  errors 
don't  return  to  zero  until  nearly  300  seconds.  04 
particular  note  is  the  constant  error  in  w  during 
saturation.  This  can  be  seen  to  be  caused  by  the  sail 
planes  since  the  error  in  w  begins  to  decrease  the  instant 
the  sail  planes  come  out  of  saturation.  Also  note  that 
although  a  25  degree  error  exists  in  pitch  (the  pitch 
command  was  -50  degrees)  at  60  seconds,  the  stern  planes  are 
only  slowing  increasing  in  deflection.  The  large  control 
surface  deflections  have  slowed  the  submarine  to  nearly  9 
knots  by  this  time,  however,  accounting  for  the  decreased 
control  gain.  Observe  that  the  sub marina  speed  rapidly 
decreases  when  the  rudder  saturates,  but  slowly  increases  in 
speed  past  100  seconds,  even  though  the  sail  planes 
saturated. 


4.10  Summary 

This  chapter  has  presented  an  overview  of  the  LQ6/LTR 
design  methodology  and  then  the  application  of  the 
methodology  to  the  design  of  a  compensator  for  the  submarine 
model . 

Compensator  designs  were  developed  for  all  four  linear 
models  of  the  submarine,  providing  the  information  necessary 


curve 


•fit 


to  apply  a  nonlinear 
parameters,  laying  the  groundwork 
algorithm  in  the  next  chapter. 


to  the  compensated 
for  the  gain  scheduling 


Time  histories  for  the  20  knot  model  were  presented  that 
further  validated  the  linearization  of  the  nonlinear 
equations  of  motion  and  displayed  the  ability  of  the 
compensator  design  to  control  the  submarine  nonlinear  model 
as  long  as  the  forward  velocity  remained  close  to  the 
nominal  point. 


Finally,  the  loss  of  control  that  can  result  from 
saturation  of  the  control  surfaces  was  demonstrated. 
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Chapter  9 


NONLINEAR  MI MO  DESIGN 


5.1  Introduction 


It  is  well  known  that  drag  and  lift  forces  on  a  body  vary 
as  the  square  of  the  velocity  1161.  Therefore,  a 
compensator  designed  for  a  particular  speed  should  not  be 
expected  to  provide  optimum  performance  over  a  widely 
varying  speed  range. 

This  chapter  demonstrates  the  technique  of  gain  scheduling 
as  a  method  of  introducing  the  inherent  nonlinearities  due 
to  velocity  into  the  compensator  design. 

The  algorithm  used  to  implement  the  gain  scheduling  is 
presented,  along  with  representati ve  nonlinear  simulations. 


5.2  The  Gain  Scheduling  Algorithm 


In  the  previous  chapters,  linear  models  of  the  submarine 
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dynamics  were  developed  for  four  different  speeds,  and  an 
LQ6/LTR  design  performed  on  each  model.  The  products  of  the 
design  effort  were  the  Kalman  Filter  and  Control  gain 
matrices  contained  in  Appendix  Cl. 

The  specifics  of  the  computer  implementation  of  the  model 
based  compensator  require  that  the  compensator  be  defined  by 
three  matrices,  A-B  G-H  C,  G,  and  H  (refer  to  equations  4.17 
and  4.18).  As  a  result,  the  gain  scheduling  algorithm  must 
produce  the  required  coefficients  to  calculate  the 
individual  entries  for  these  matrices. 

The  fact  that  lift  and  drag  forces  vary  as  the  square  of 
the  velocity  indicates  the  choice  of  a  polynomial  least 
squares  algorithm  to  produce  the  required  coefficients. 
Additionally,  application  of  a  least  squares  algorithm  to 
matrices  as  large  as  11  x  11  requires  that  a  computer  be 
used.  The  fortran  program  in  Appendix  D1  is  an  adaptation 
of  the  vector  based  program  developed  by  A.  Miller  C17J  that 
will  accommodate  square  and  non-square  matrices  up  to  order 
11. 

The  program  calculates  correlation  coefficients  for  each 
coefficient  produced,  providing  a  determination  of  how  good 
the  appr ox i mat i on  is.  As  expected,  a  linear  fit  (polynomial 
of  order  1)  produced  very  poor  correlation.  A  second  order 
curve  fit  produced  correlation  coefficients  averaging  B.98, 


with  a  third  order  -fit  only  marginally  better.  Since 
implementation  of  gain  scheduling  would  require  187  *  n 
multiplications  and  187  *  n  additions,  where  n  =  the  order 
of  the  polynomial,  a  minimal  order  polynomial  is  desired. 
Thus  a  second  order  polynomial  was  used  to  calculate  the 
coefficients.  The  results,  along  with  the  correlation 
coefficients,  are  presented  in  Appendix  D2. 


The  fortran  program  implementing  the  submarine  model  at 

CSDL  was  then  modified  to  permit  the  coefficients  for  the 

A-B  G-H  C,  6,  and  H  matrices  to  be  read  in,  and  new  values 

for  these  matrices  calculated  for  each  iteration  during 

integration  of  either  the  linear  or  nonlinear  equations  of 

motion.  To  calculate  a  particular  entry  for  one  of  the 

matrices,  say  the  '.i,j)  element  of  the  H  matrix,  the  current 

forward  velocity  u(t)  is  first  converted  to  knots  <u^>. 

Then,  if  c.  .  c.  .  and  c.  .  _  are  the  coefficients 
i,J,u>  i  ,J*i  x,j,^ 

calculated  by  the  gain  scheduling  algorithm  for  that  element 
of  the  H  matrix,  the  updated  value  is  given  by 


*  c.  .  _  +  u.  <c.  .  ,  +  u.  (c.  .  _>  > 

l , j,0  k  i  ,  j  ,  1  k  l  ,  j  ,2 


(5.1) 


Note  that  the  order  of  multiplication  and  addition  has  been 


arranged  to  require  the  minimum  amount  of  calculation. 


The  fact  that  the  computer  implementation  of  the  submarine 
model  did  not  provide  for  propeller  rpm  to  be  varied  was 
overcome  by  setting  the  rpm  to  the  final  speed  and 


introducing  the  initial  speed  into  the  initial  conditions 


5.5  Testing  of  the  Nonlinear  Compensator 

Although  the  above  method  of  overcoming  the  constant  rpm 
constraint  of  the  model  permitted  simulation  of  gain 
scheduling,  the  flexibility  of  the  simulation  was  somewhat 
limited.  Only  one  speed  excursion  could  be  accomplished 
during  any  single  run,  with  the  propeller  always  rotating  at 
a  constant  rpm. 

Figure  5.1  simulates  changing  depth  without  imparting  a 
pitch  angle  to  the  submarine  and  provides  an  excellent 
example  of  both  the  usefulness  of  w(t>  as  a  state  variable 
and  the  increase  in  control  surface  gain  produced  by  the 
gain  scheduling  algorithm.  The  time  history  of  the  sail 
planes  (db)  indicates  their  dominant  role  in  this  maneuver. 
Although  the  command  in  w(0.5  ft/sec)  is  attained  quite 
early  at  60  seconds,  the  submarine  is  still  slowing  and  the 
nonlinear  compensator  increases  the  sail  plane  deflection 
accordingly.  The  second  command  in  w,  at  275  seconds, 
produces  the  same  magnitude  error  signal,  but  results  in  a 
much  larger  initial  deflection  of  the  sail  planes  due  to  the 
lower  speed.  The  error  in  w  is  seen  to  stay  constant  and 
very  near  zero  during  the  sail  plane  movement. 
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<f)  Fair water  plans  deflection 


(g)  Rudder  deflection 


Figure  5.1.  (Cont'd) 
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The  next  two  simulations  reflect  an  attempt  to  further 
define  the  nature  of  the  cross  coupling  that  exists  between 
pitch  angle  0  and  yaw  rate  r.  Both  tests  have  initial 
conditions  corresponding  to  the  nominal  point  for  the  5  knot 
model  with  the  final  speed  equal  to  the  nominal  speed  for 
the  25  knot  model.  Command  inputs  were  applied  as  step 
signals  that  corresponded  to  the  nominal  values  for  the 
three  output  variables  for  the  10,  20,  and  25  knot  models. 
The  step  inputs  for  the  10  and  20  knot  nominal  points  were 
timed  so  as  to  occur  simultaneously  with  the  submarine 
attaining  the  forward  velocity  appropriate  to  that  nominal 
point. 

In  figure  5.2,  the  submarine  is  allowed  to  settle  for  a 
short  time  at  the  20  knot  nominal  conditions  before  the 
final  step  change  is  applied.  In  figure  5.3,  however,  the 
last  step  change  is  applied  as  soon  as  the  appropriate 
velocity  is  reached.  Thus  the  two  time  histories  are 
identical  up  to  about  35  seconds.  A  comparison  of  the  two 
sets  of  time  histories  reveals  that  the  cross  coupling 
between  yaw  rate  and  pitch  angle  is  not  minimized  in  the 
second  test,  although  this  might  be  expected  since  the 
submarine  should  be  experiencing  an  attitude  very  similar  to 
that  for  which  the  compensators  were  designed.  Instead, 
although  the  curves  appear  smoother  in  the  second  test,  a 
close  inspection  of  the  errors  and  the  control  surface 
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(•)  Stern  plane  deflection 


(f)  Fairwater  plane  deflection 
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deflections  reveal  that  the  degree  o-f  cross  coupling  is  a 
•function  o-f  the  yaw  rate  error  at  the  time  when  a  step 
change  is  applied.  In  figure  5.2(c),  the  error  in  the  yaw 
rate  when  the  third  step  is  applied  is  approximately  (3.06)3 
rad/sec  whereas,  in  figure  5.2(c),  it  is  nearly  twice  as 
large.  The  result  is  a  larger  overshoot  for  all  three 
outputs  and  an  increased  rudder  deflection.  Furthermore, 
although  the  third  step  change  was  applied  50  seconds 
earlier  in  the  second  test,  the  errors  at  200  seconds  are 
larger  than  for  the  first  test. 

Time  did  not  permit  further  investigation  of  the  cross 
coupling  present  in  the  submarine  model  simulations.  It 
should  prove  informative,  however,  to  compare  the 
simulations  of  figures  5.1,  5.2,  and  5.3  to  identical 
simulations  for  submarine  models  linearized  about  a  straight 
and  level  trajectory. 

The  last  test,  shown  in  Figure  5.4,  illustrates  a  speed 
increase  from  5  to  25  knots  with  all  commands  set  to  zero. 
This  simulation  reveals  an  oscillatory  response  not  present 
in  any  of  the  simulations  presented  so  far.  It  occurred 
consistently  for  all  tests  of  the  nonlinear  compensator 
involving  a  speed  increase  to  approximately  30  knots  and 
having  significant  time  for  settling  out  to  steady  state. 
It  should  be  noted  that  the  oscillations  are  quite  small  in 


amplitude  and  are  very  lightly  damped. 


Although  the  cause  of  the  oscillations  was  not  determined , 
some  phase  lag  appears  to  be  present  in  the  stern  plane 
deflection  compared  to  the  error  in  6.  This  indicates  the 
requirement  for  more  data  points  to  accurately  predict  the 
model's  behavior  at  higher  speeds. 


5.4  Summary 

This  chapter  discussed  the  motivation  for  including  the 
nonlinear  aspects  of  the  submarine  model's  response  into  the 
compensator  design.  The  technique  of  gain  scheduling  was 
introduced  as  the  method  by  which  this  could  be 
accomplished,  and  the  specific  algorithm  used  was 
presented. 

The  first  gain  scheduling  simulation,  figure  5.1, 
highlights  the  application  of  this  technique  in  the  time 
history  of  the  sail  planes  (figure  5.1(f)).  Simulations 
using  the  full  nonlinear  controller  were  also  presented  that 
further  demonstrated  the  cross  coupling  present  in  the 


submarine  model 


t-'iJ 


(a)  St arn  plans  dafl action 


VA 


2*0  300 

TIME 


< •f >  Fairwatar  plana  da-fl action 


240  303 

TIME 


(g)  Ruddar  da-flaction 


r-r^'t  ^“#.*r**  .  ■»*  • .  v. 


Figura  5.4.  (Cont'd) 
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Chapter  b 


SUMMARY  AND  DIRECTIONS  FOR  FURTHER  RESEARCH 


6. 1  Summary 


This  thesis  has  presented  a  mul ti variable  control  design 
example  consisting  of  the  following  steps: 

1.  Four  linear  models  of  a  full  size  submarine  were 
developed  using  a  computer  implementation  of  the 
linear  and  nonlinear  equations  of  motion. 

2.  The  resulting  models  were  reduced  in  order  by 
inspection  of  the  A  and  B  matrices  and  subjected  to 
modal  analysis  to  determine  their  suitability  for  the 
LQG/LTR  design  procedure. 

3.  Response  characteristics  of  the  open  loop  nonlinear 
model  were  used  to  establish  guidelines  for 
performance  specifications  and  additional  dynamics 
were  included  in  the  models  in  the  form  of  integrators 
in  each  input  channel. 


4.  The  controller  was  designed  using  the  LQG/LTR  combined 
time  and  frequency  domain  methodology  as  formulated  by 
Stein,  Doyle,  et.al.  Loop  shaping  techniques  were 
independently  developed  that  provided  identical 
behavior  of  the  singular  values  at  both  ends  of  the 
frequency  spectrum  for  linear  models  using  integral 
control . 

5.  A  nonlinear  capability  was  added  to  the  controller 
design  using  the  technique  of  gain  scheduling.  A 
second  order  polynomial  was  used  to  produce 
coefficients  that  enabled  prediction  of  the 
compensator  matrices  for  nominal  points  not  modeled. 


6.2  Some  Conclusions  and  Directions  for  Further 
Research 


Modal  analysis  should  be  used  as  a  first  step  in  any 
mul ti variable  controller  design.  It  not  only  provides  all 
of  the  necessary  information  for  proceeding  with  the  design, 
but  also  provides  invaluable  reference  material  for 
determining  the  cause  of  any  anomalies  in  the  behavior  of 
the  model . 


The  techniques  used  in  this  thesis  to  match  the  singular 


values  of  the  Kalman  Filter  transfer  matrix  are  simple  and 
straight-f or ward.  Unfortunately,  they  only  apply  when 
certain  of  the  matrices  used  are  invertible.  An  extension 
of  this  method  to  cover  the  case  of  noninvertible  matrices 
using  a  perturbation  matrix  would  provide  a  useful  and 
general  tool  to  the  control  system  designer. 

The  extreme  usefulness  of  singular  value  characterizations 
for  multivariable  systems  was  demonstrated  in  this  thesis. 
The  concept  is  an  easily  understood  extension  of  the 
familiar  Bode  plots  and  provides  a  convenient  means  of 
describing  and  then  ensuring  the  attainment  of  performance 
requirements  for  mul ti variable  systems. 

The  cross  coupling  between  yaw  rate  and  pitch  angle 
displayed  by  the  submarine  models  used  in  this  thesis 
requires  further  investigation.  In  particular,  performance 
and  transient  response  for  models  linearized  about  a 
straight  and  level  trajectory  should  be  compared  to  the 
simulations  presented  herein.  Such  a  comparison  should 
allow  a  determination  to  be  made  as  to  whether  or  not  gain 
scheduling  needs  to  performed  on  the  basis  of  roll  angle  as 
well  as  forward  velocity. 

The  gain  scheduling  algorithm  described  in  the  previous 
chapter  is  an  idealization  in  that  the  compensator  matrices 
were  updated  on  each  iteration  of  the  integration.  An 


actual  implementation  of  gain  scheduling  might  require  that 
an  overlapping  step  type  of  approach  be  used  to  decrease 
both  the  computational  burden  and  the  possibility  of 
introducing  instabilities  into  the  controller. 

In  so  far  as  the  purpose  of  this  thesis  was  to  present  an 
LQG/LTR  design  example  and  demonstrate  the  technique  of  gain 
scheduling,  the  reference  plane  problem  was  ignored. 
Unfortunately,  controlling  a  submarine  in  an  inertial 
reference  frame  utilizing  a  body  coordinate  system  somewhat 
abstracts  the  state  space  description  of  the  submarine 
model.  To  make  the  conversion  externally,  however,  presents 
the  control  system  designer  with  a  time  varying  C  matrix. 
Replacement  of  the  state  variables  w  and  r  with  the 
derivatives  of  z  and  ¥  would  provide  a  more  useful  set  of 
state  variables  for  the  purposes  of  controller  design. 

Although  the  submarine  has  six  degrees  of  freedom  in  its 
environment,  the  three  control  surfaces  permit  the  control 
of  only  three  output  variables.  Of  course  the  constant  rpm 
constraint  placed  on  the  model  in  this  thesis  is  artificial, 
and  was  intended  only  to  reflect  current  operating 
procedures.  The  addition  of  propulsion  control  and 
differential  action  for  the  control  surfaces  (rudders  and 
stern  planes)  would  provide  a  total  of  6  controls,  yielding 
a  much  more  flexible  system. 


Appendix  A 


State  Space  Matrices  for  tha  Linaarizad  Modal a 


The  elements  of  the  A  and  B  matrices  are  presented  in  the 
standard  row  and  column  format.  In  the  case  of  the  A 
matrices,  the  10  elements  of  each  row  are  displayed  as  two 
rows,  containing  the  first  five  elements  in  one  row  and  the 
second  five  elements  in  the  next. 


Model  S5R5 


A  Matrix 


-.8532E-02 
.  1590E+01 

. 1532E-02 
2569E+01 

2690E-03 

1840E-01 

-.3420E-04 

1544E-01 

- . 7759E— 05 
-.5396E-03 

2941E-04 

8534E-01 

. 0000E+00 
2027E— 01 

. 0000E+00 
. 1050E-01 

. 0000E+00 
. 1000E+01 

. 2027E— 01 
. 0000E+00 


1 166E-01 
. 0000E+00 

4036E-01 
.  1318E+00 

. 6475E-03 
. 0000E+00 

2184E-02 

1621E+00 

. A245E-05 
. 0000E+00 

2764E-03 
. 2618E-03 

. 0000E+00 
.5899E-12 

. 0000E+00 

. 6547E-02 

. 0000E+00 
291 IE— 10 

1050E-01 
. 7882E+00 


1122E-02 

. 0000E+00 
-.7633E-01 
. 0000E+00 
1845E-02 
. 0000E+00 
. 0000E+00 
. 0000E+00 


1594E-03 
. 2954E-03 

. 1187E-02 
2805E-04 

2266E-01 
. 7616E-02 

7419E-03 
. 3450E-04 

. 2208E-03 
2512E-02 

7509E-05 
— . S572E— 07 

. 0000E+00 
6549E-02 

. 0000E+00 
. 0000E+00 

. 0000E+00 

. 1327E-03 

. 9997E+00 
— .  9643E+-01 

B  Matrix 

— . 2826E-02 

. 0000E+00 
5&24E-01 
. 0000E+00 
. 5402E-03 
. 0000E+00 
. 0000E+00 
. 0000E+00 


1001E-01 
. 0000E+00 

3551E+00 
. 0000 E+ 00 

8010E+00 
. 0000E+00 

8579E— 01 
. 0000E+00 

-.6028E-02 
. 0000E+00 

134&E-02 

. 0000E+00 

. 1000E+01 
. 0000E+00 

. 0000E+00 
. 0000E+00 

. 0000E+00 
. 0000E+00 

. 0000E+00 
. 0000E+00 


U37E-01 
.8949E-01 
2885E-06 
. 1661E-02 
. 9514E-07 
-.226&E-02 
. 0000E+00 
. 0000E+00 


. 203AE+00 
.  0000E+00 

.  1 268E+00 
.  0000E+00 

.  2183E+01 
. 0000E+00 

1323E+00 
. 0000E+00 

- .  8562E-0 1 
. 0000E+00 

—.71 10E— 04 
. 0000E+00 

.2128E-03 

. 0000E+00 

. 9999E+00 
. 0000E+00 

1050E-01 
. 0000E+00 

. 0000E+00 
. 0000E+00 


Model  S10R5 


A  Matrix 


1491E-01 

20&7E-01 

7935E-03 

1777E-01 

2832E+0 1 

. 0000E+00 

.  2948E-03 

.  0000E+00 

2771E-02 

7168E-01 

. 1182E-02 

— .  5860E+00 

45&2E+01 

. 1314E+00 

3061E-03 

.  0000E+00 

4935E-03 

. 1470E-02 

— .  4014E-01 

1427E+01 

1274E-01 

. 0000E+00 

.  7599E-02 

.  0000E+00 

6136E-04 

— . 3913E-02 

1 330E-02 

1522E+00 

2096E-0 1 

1616E+00 

. 3765E-03 

. 0000E+00 

1467E-04 

. 6635E-05 

. 3919E-03 

1069E-01 

1 162E-02 

. 0000E+00 

2506E-02 

. 0000E+00 

5373E-04 

— . 4895E-03 

-.7697E-05 

— . 2670E— 02 

1515E+00 

•  26UE-03 

— . 608 IE— 06 

. 0000E+00 

0000E+00 

. 0000E+00 

. 0000E+00 

. 1000E+01 

7024E-01 

. 981 BE- 10 

U68E-01 

. 0000E+00 

0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

3312E-01 

.  1162E-01 

.  0000E+00 

. 0000E+00 

0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

1002E+01 

1400E-08 

.  8189E-03 

. 0000E+00 

7011E-01 

— .  3304E-01 

.  9970E+00 

. 0000E+00 

0000E+00 

.  1398E+01 

1709E+02 

. 0000E+00 

B  Matrix 

— .  3512E— 02 

8903E-02 

— . 3582E— 01 

. 0000E+00 

e  0000E+00 

.  2910E+00 

2467E+00 

1772E+00 

— .  9087E— 06 

.  0000E+00 

. 0000E+00 

.  5191E-02 

5992E-02 

.  1702E-02 

. 2997E-06 

. 0000E+00 

.  0000E+00 

— .  7404E-02 

.  0000E+00 

.  0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

a  0000E+00 

. 0000E+00 

. 2705E+00 
. 0000E+00 

. 2069E+00 
. 0000E+00 

.3871E+01 
. 0000E+00 

— . 2374E+00 
.  0000E+00 

151 9E+00 
. 0000E+00 

. 8318E-03 
. 0000E+00 

. 232BE-02 

. 0000E+00 

.  9995E+00 
. 0000E+00 

3320E-01 

.  0000E+00 

. 0000E+00 
. 0000E+00 


Model  S20R5 


A  Matrix 

-.2889E-01 
. 5643E+0J 

-.4046E-01 
. 0000E+00 

-.&327E-02 
. 2754E-03 

3487E-01 
. 0000E+00 

. 5539E-02 
9060E+01 

1416E+00 
.  1216E+00 

8891E-02 

-.6676E-02 

-.6934E-M80 
. 0000E+00 

1206E-02 
. 1720E+00 

. 7565E-02 
. 0000E+00 

8581E-01 
. 7099E-02 

-.2851E+01 
. 0000E+00 

— . &092E-04 
.2341E-01 

8099E-02 

1496E+00 

— . 2414E-02 
.821  IE-02 

- . 3005E+00 
. 0000E+00 

-.3405E-04 

5757E-02 

— . 3369E-04 
. 0000E+00 

• 7281E-03 
— . 2342E-02 

-.2096E-01 
. 0000E+00 

1078E-03 

2994E+00 

--9A20E-03 
- 2417E-03 

. 4209E-04 
1326E-04 

-.792BE-02 
. 0000E+00 

. 0000E+00 
— . 3853E+00 

. 0000E+00 
9059E-09 

.  0000E+00 
—•  2650E-01 

. 1000E+01 
.  0000E+00 

. 0000E+00 
. 1 397E+00 

•  0000E-f-00 

. 2302E-01 

. 0000E+00 
. 0000E+00 

.  0000E+00 
.  0000E+00 

. 0000E+00 
.  1063E+01 

. 0000E+00 
. 2498E-08 

. 0000E+00 
. 9609E-02 

. 0000E+00 
. 0000E+00 

.  3&2&E+00 
. 0000E+00 

1 302E+00 
• 2619E+01 

. 9228E+00 
— . 3159E+02 

. 0000E+00 
. 0000E+00 

B  Matrix 

1372E-01 

. 0000E+00 
9636E+00 
. 0000E+00 
2341E-01 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+00 

-.347BE-01 
. 0000E+00 
-.6921E+00 
. 0000E+00 
. 6648E-02 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+00 

1399E+00 
. 1137E+01 
-.3550E-05 
. 2028E-01 
. 1171E-05 
-.2892E-01 
. 0000E+00 
• 0000E+00 
. 0000E+00 
. 0000E+00 

— .  4169E+-00 
.  0000E+00 

. 1803E+00 
. 0000E+00 

.7576E+01 
. 0000E+00 

-.4355E+00 
. 0000E+00 

-.299BE+00 
.  0000E+00 

. 1271E-01 
.  0000E+00 

• 5434E-0 1 
. 0000E+00 

.  99 02E+00 
. 0000E+00 

1499E+00 
. 0000E+00 

. 0000E+00 
. 0000E+00 


Model  S25R5 


A  Matrix 


3521E-01 

4803E-01 

1255E-01 

-.4146E-01 

1606E+01 

. 6762E+01 

. 0000E+00 

. 2343E-03 

. 0000E+00 

. 0000E+00 

. &953E-02 

1718E+00 

2605E-01 

3137E+00 

1038E+00 

1124E+02 

. 1017E+00 

1850E-01 

. 0000E+00 

. 0000E+00 

1901E-02 

. 1686E-01 

1 060E+00 

-.3427E+01 

. 9191E+01 

. 4173E+00 

. 0000E+00 

. 6039E-02 

. 0000E+00 

. 0000E+00 

. 3581E-04 

1038E-01 

2127E-02 

- . 3679E+00 

3943E+00 

. 1083E+00 

1 251E-*-00 

• 2275E-01 

. 0000E+00 

. 0000E+00 

-.4478E-04 

— . 9358E— 04 

. 8713E-03 

-.2491E-01 

— .  368 IE +00 

132SE-01 

. 0000E+00 

1992E-02 

. 0000E+00 

.  0000E+00 

1331E-03 

1 197E-02 

. 1 125E-03 

1242E-01 

•  29B5E-01 

-.3652E+00 

.  2020E-03 

-.367&E-04 

. 0000E+00 

• 0000E+00 

. 0000E+00 

.  0000E+00 

.  0000E+00 

. 1000E+01 

. 1770E+00 

-.7482E+00 

2864E-08 

-.4430E-01 

.  0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

.  0000E+00 

.  0000E+00 

. 9731E+00 

. 2303E+00 

. 2784E-01 

. 0000E+00 

.  0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

2905E+00 

. 1227E+01 

. 4699E-08 

. 2700E-01 

. 0000E+00 

. 0000E+00 

. 6095E+00 

1 825E+00 

.7715E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 2750E+01 

3285E+02 

. 0000E+00 

. 0000E+00 

B  Matrix 

— • 2061E— 01 

— • 5221E-01 

2101E+00 

.0000E+00 

. 0000E+00 

. 1703E+01 

1444E+01 

1039E+01 

5329E-05 

. 0000E+00 

. 0000E+00 

. 3046E-01 

3508E-01 

. 9981E-02 

. 1 758E-05 

.  0000E+00 

. 0000E+00 

4333E-01 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

. 0000E+00 

.  0000E+00 

. 0000E+00 

. 0000E+00 
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All  matrices  are  presented  in  the  standard  row  and  column 
-format.  Additionally,  the  data  in  appendices  B2  through  B6 
are  complex  numbers.  These  numbers  are  always  displayed 
with  the  imaginary  part  directly  below  the  real  part.  For 
example,  each  row  of  a  matrix  with  complex  entries  would  be 
displayed  as  two  rows,  the  real  parts  in  one  row  and  the 
corresponding  imaginary  parts  in  the  second  row.  The 
eigenvectors  are  presented  as  (complex)  column  vectors. 


Appendix  81 


Matrices  Used  to  Pertore  Unit  Transtoreations 


1.  Matrix  used  to  preaultiply  both 

the  A  and  8 

eat n cast 

l.aaiaEsN 

a.aaaae^H 

i.mhe+ii 

e.eeaeE+ee 

a.eeeeE+ea 

B.NNE^N 

e.eeeeE+ee 

a.aaeeE+ea 

l.eaaaE+ee 

a.aeeiE+aa 

e.eaeeE+ea 

8. 8»»8€>I0 

0.BBBBE+0B 

a. aaatE^aa 

a.eeeaE+ea 

e.eaeeE+ea 

l.eeeeE+ea 

e.aeaaE+ae 

e.eeeaE+ee 

a. aaaaE^aa 

e.eaeeE+ea 

a.eeaaE+aa 

0.00ME+00 

e.eeaeE+ee 

5.7388E+01 

e.eeeaE+ee 

e.eeeeE+ee 

e.eeeeE+ee 

a.eeeaE+aa 

e.eaeeE+ea 

e.eeaeE+ee 

a.eeeeE+ea 

5. 73B8E+-01 

e.eeeeE+ee 

e.eeeeE+ee 

0. 0000E+0B 

0.0000E+00 

0.000IE+00 

0.0000E+00 

0.0000E+00 

5.730IE+01 

e.eeeeE+ee 

0.0000E+00 

0.0000E+00 

0.0000E+00 

0.0000E+00 

0.0000E+00 

0.0000E+00 

5.7300E+01 

0.0000E+00 

0.0000E+00 

0.0000E+00 

0.0000E+00 

0.8808E+00 

0.0000E+00 

0.0080E+00 

!.  Matrix  used  to  postauitiply  the 

A  aatrix  only: 

l.aeaaE+ea 

B.0000E+0B 

B.0000E+00 

8.00B0E*08 

e.eieeE+ea 

0.0000E+00 

0.0000E+08 

a.BaaeE+aa 

1.0000E+00 

e.aeaaE+ae 

a.aeaaE+aa 

a.aeeaE-sfle 

e.eeeeE+ee 

0.0080E+80 

0.0000E+00 

0.0000E+00 

1 . 0000E+00 

0.0000E+00 

0.0000E+00 

0.0000E+00 

0.0080E+80 

e.aaaaE+aa 

a.aaaaE+ea 

8.8000E+08 

1 . 7452E-02 

B.B000E+00 

0.0000E+00 

e.eeeeE+ee 

a.aaaaE+aa 

0.0000E+00 

e.eaeeE+aa 

a.eeeeE+ea 

1.7452E-02 

e.eeeeE+ee 

0.000iEt0i 

a.aaaaE+ea 

0.0000E+08 

e.eeeaE+ee 

e. B8B8E*aa 

a.eeeeE+ea 

1 . 7452E-B2 

e.eeeeE+ee 

a.aaaaE+ea 

0.0888E*00 

aitaiEM 

8i}80EM 

e.eeeaE+ee 

e.eeeeE+ee 

1 . 7452E-B2 

i.eauE+ee 

e.eeeaE+ee 

a.aaaaE+ea 

l.88ME+8a 

a.eeeeE+ea 

e.eeeeE+ee 

8.800K*0B 

3.  Matrix  used  to  postauitiply  the  8  aatnx: 

S.73N£s|i  8.0000E+00  0.0000E+00 

I.MM£*M  5.73ME+II  I.NNE^H 
I.NNE^N  I.NNEsN  5.7300E*B1 


I.  MIIE+M 
0.0B0BE+00 
a.aaaaE+ea 
a.aaaaE+aa 
e.eaaaE>ea 
a.aaaaE+ea 
a.eaeeE^ea 

5.7300E+81 

a.eaaaE^aa 

a.aaaeE+aa 

0.000iE*aa 

0.0000E+00 

i.eeeeE+ea 

a.aaeaE+ea 

0.0000E*00 

1.7452E-B2 


Appendix  B2 


*1. 1277E-02 

0.0000E*00 


9.9281E-01 

0. 0000E+00 

1. 1961E-01 

0.BB00E+00 

■1.8966E-03 

B.0000E+00 

-2.9432E-05 

0.0000E1-00 

7.o272E-05 

8.0000E*00 

-8.4622E-04 

0.00006*00 

-1.7819E-03 

0.0000E+00 

•4.9409E-03 

0.0000E00 


Modal  Decoeoosi t i on  far  Model  S5R5 


Eigenvalues 


-2.4709E-02  -2.7689E-02 
0.I000E+00  1.7643E-02 


7.4170E-01  8.9904E-02 
0. 8000E+00  3.2552E-02 

6.6500E-01  2. 9770E-02 
0. 0000E+00  -7. 3528E-03 

B.59E2E-02  7.5459E-01 

0.80886*00  5.9693E-01 

2.5411E-04  1.7394E-03 
0.0880E*00  1.3394E-04 

-2.4696E-04  -B.8074E-03 
0.0000E+00  4.8192E-03 

-3.5640E-83  -3. 1491E-04 
0.0000E+00  -2.7189E-06 

-9.0&14E-03  1.1554E-03 
0.0000Et00  -6.8U6E-03 

1.2868E-02  2.5288E-01 
0.0000E+00  -1. 1295E-02 


-2.7689E-82  -4. 1686E-82 
-1.7643E-02  3.97ME-01 


Eigenvectors 

8.O904E-02  -2.9654E-03 
-3.2552E-02  -3.2693E-03 

2.9770E-02  -2.20S6E-01 
7.3528E-03  -2.7560E-01 

7.5459E-01  -5.5B26E-01 
-5.9693E-01  -5. 4402E-02 

1.7394E-03  -5.1519E-02 
-1.3394E-04  2.7309E-01 

-6.8074E-03  -4.0647E-03 
-4.8192E-03  -9.2717E-04 

-3. 1491E-04  -5.9911E-04 
2.7189E-06  -S.6137E-84 

1. 1554E-03  6. 9284E-01 
6.8116E-03  5.687BE-02 

2.5288E-01  -1.5238E-03 
1. 1295E-02  -I.0097E-03 


-4. 1688E-82  -5.3308E-02 
-3.9761E-01  0.0808E*B8 


-2.9i5iE-03  -e.5228E-02 
3.2693E-03  0.0000E+00 

-2.20B6E-01  1.M24E-02 
2.75602-01  0.0000E+00 

-5.5B26E-01  -9.6303E-01 
5. 4402E-02  0.0000E+00 

-5.1519E-02  -1.3353E-03 
-2.7309E-01  0.0000E*00 

-4.0647E-03  1.3928E-02 
9.2717E-04  0.0000E*00 

-5.9911E-04  1.1526E-04 
S.6137E-04  0.0000E+00 

6.9284E-01  -6.9590E-03 
-5.6878E-02  0.0080E+00 

-1.5238E-03  -2.6043E-01 

1.0097E-03  0.0000E+00 


-9.8265E-02 

0.0000E+00 


2.4125E-01 

0.0000E+00 

-9.7007E-01 

0.0000E*00 

-B.9253E-03 

0.0000E+00 

-1.B784E-03 

0.0000E*00 

-2.5031E-04 

0.0000E*0i 

-2.0701E-02 

0.0000E+00 

1.5096E-02 

8. 0810E*00 

3.7533E-03 

0.0000E+00 


$ 


Appendix  62  (cont  d) 
Model  S5R5 


Controllability  Matrix 

3.2358E-83  2.2815E-83  -2.3825E-01 
1.8821E-18  -1 . 2869E-10  1.358BE-11 

8.3133E-83  -8.5145E-83  2.8357E-81 
4.5845E-11  6. 1762E-I l  -1.375BE-11 


-8.4924E-82  5.7778E-82  -7. 971  IE-03 
8.5281E-82  2.7335E-82  5.7777E-83 

-8.4924E-82  5.7778E-82  -7. 971  IE-83 
-8.5281E-82  -2. 7335E-82  -5.7777E-83 


-1.4243E-83  1.5582E-84  4.984BE-B4 
7.2981E-85  4.1261E-05  -3.295BE-B3 


-1.4243E-B3  1.5582E-B4  4.9B4BE-B4 
-7.29B2E-B5  -4. 1261E-B5  3.295BE-83 

-1.5716E-81  1. 141BE-B1  4.5351E-B3 
-7.2441E-18  -6.8431E-1B  -1.7913E-1B 

2.5416E-B4  2.5224E-B4  7.B342E-B2 
7.8B13E-12  7.6491E-13  3.6I12E-13 


Observability  Matrix 

1 . 89AAE-B4  B.5942E-B3  7.S459E-B2  7.5459E-B2  -5.5826E-I2  -5.5824E-B2  -9.43B3E-B2  -B.9253E-B4 

I.BBHE^BB  B.8B8BE+88  5.9493E-B2  -5.9493E-B2  -5.44S2E-B3  5.44B2E-B3  B.BIBBE+BB  B.BBIBE*B8 

B.4422E-B4  -3.544BE-B3  -3.1491E-B4  -3.1491E-B4  -5.9911E-B4  -5.99UE-B4  1.1524E-I4  -2.B7B1E-B2 

B.SBBBE+BB  B.SBBBE+H  -2.7189E-B4  2.71B9E-B4  -8.6137E-I4  8.4137E-B4  fl.BBSBE+Bfl  1. 08I0E+8B 

4.94B9E-B3  1.2848E-B2  2.S288E-B1  2.5288E-B1  -1.5238E-B3  -1.5238E-B3  -2.4B43E-B1  3.7533E-B3 

B.BME+BB  I.IBBBE*H  -1.1295E-B2  1.1295E-B2  -1.N97E-I3  1.H97E-03  0.B0I8E+I8  B.BIBBE+BB 


Appendix  83 


-1.9940E-02 

0.00008+00 


-9.9044E-01 

0.0000E+00 

-1.2375E-01 

0.0000E+00 

-5.9146E-02 

0. 00008*00 

1.437SE-04 

0.0000E+00 

-3.5989E-04 

0.00008*00 

B. 71448-04 
0. 0000E+00 

4. 1896E-03 

0. 00008+00 

1.4149E-02 

1.00008+00 


Modal  Decomposition  for  Model  S10R5 


Eigenvalues 


-2.2087E-02  -2.2087E-02 
1.5563E-02  -1.5563E-02 


■4.7861E-02  -7.2735E-02  -7.2735E-02 
0.0000E+00  3.B643E-01  -3.8643E-01 


■1.4962E-01  -1.7528E-01 

0.80008+00  0.00008+00 


Eigenvectors 


-5.4'e2E-03 

-5.4762E-03 

-7.3665E-01 

-9.3156E-03  -9.3156E-03 

6.23798-02 

2.3753E-31 

4.9564E-02 

-4. 9564E-02 

0.00808+00 

4.3473E-03  -4.8473E-03 

0.0000E+00 

0.0000E+00 

2.8387E-03 

2.0387E-03 

-6.7361E-01 

-2.9934E-01  -2.9904E-0! 

-3.5011E-02 

-9.7065E-01 

4.3490E-32 

-4.3490E-02 

0.0000E+00 

2.2I70E-01  -2.2070E-01 

8.80008+80 

8.00008+00 

6.5223E-01 

6.5223E-01 

-5.6882E-02 

-2.0650E-B!  -2.08508-01 

9.7768E-01 

4.7148E-04 

7.2693E-01 

-7.2698E-01 

0.0000E+00 

7.1837E-01  -7. 1837E-01 

0.0600E+00 

0.0000E+00 

1.5976E-03 

1.5976E-03 

-6. 46658-04 

1.8240E-01  1.8240E-01 

-3. 1077E-03 

-6. 2008E-03 

-1.S416E-83 

1.S616E-03 

0.00008+00 

8.5333E-02  -8.5333E-02 

8.30008+00 

0.0800E+00 

5.012AE-04 

5.0126E-04 

6.4362E-05 

-2.61B0E-03  -2.61B0E-03 

-2.9157E-02 

-1.9973E-03 

5.6S92E-03 

-5.6592E-03 

0.0000E+08 

4.7088E-03  -4.7088E-03 

0.0000E+00 

0.0000E+00 

-1.0958E-04 

-1.0958E-04 

3.6271E-03 

-1.3168E-03  -1.3168E-03 

-2.0260E-04 

-2.0382E-02 

-1.6359E-04 

1 . 6859E-04 

0.00008+00 

5.2738E-04  -5.2738E-04 

0.0000E+00 

0.0000E+00 

-7.3594E-03 

-7.3594E-03 

1 . 6890E-02 

1.2725E-01  1.2725E-01 

3.6U7E-B2 

2.8127E-02 

-1.4405E-02 

1.440SE-02 

0.00008+00 

-4.9625E-01  4.9625E-01 

0.0000E+00 

0.0000E+00 

1 . 0443E-01 

1.0443E-01 

-7.9559E-03 

-2. 0542E-03  -2.0542E-03 

1.92B0E-01 

1 . 3375E-02 

-1.7467E-01 

1.7467E-01 

0.0000E+00 

3.4431E-03  -3.4431E-03 

B.0000E+00 

0.0000E+00 

1 

1 


Appendix  B3  icont'd! 


Model  S10R5 


Controllability  Matrix 

2.3370E-02  2. 9A72E-83  7.7715E-I1 
-8. 5265E-1 1  1.2147E-1I  -1.1231E-10 

-2. &658E-81  -1.5415E-02  -2.3874E-02 
4.0744E-02  5.3014E-02  -3.6798E-02 

-2.6650E-81  -1.5415E-82  -2.3874E-02 
-6.8764E-02  -5. 301 4E-82  3.6790E-02 

-1.5316E-02  -5.2494E-03  -9.1870E-01 
-1 . 4236E-1B  8.9197E-11  5.8874E-11 

-2.1579E-83  -4.6727E-84  1.7865E-02 
•1. 1457E-02  1.304 IE-03  2.4127E-03 

-2. 1579E-03  -4.6727E-84  1.7865E-02 
1.1457E-02  -1.3061E-03  -2.4127E-03 

1 . 7643E-01  -7.9939E-02  -1.B003E-82 
-9.0614E-U  3.3785E-10  -1.8444E-10 

1.2407E-03  1.0258E-03  2.3037E-01 
-9.7488E-12  1.0880E-11  1.4121E-13 


Observability  Matrix 

-5.9144E-I3  4.5223E-02  6.5223E-02  -5.4862E-03  -2.0650E-02  -2.0650E-02  9.7748E-02  4.7148E-05 

0.0000E+00  7.2698E-02  -7.2498E-02  0.8000E+00  7.1337E-02  -7.1837E-02  0.0000E»00  8.80888*80 

B.7144E-04  -1.095BE-04  -1.0958E-04  3.4271E-03  -1.3168E-03  -1.3168E-03  -2.0280E-04  -2.0382E-02 

0. §88>E*00  -1.4859E-04  1.4859E-04  8. 8000E*I0  5.2738E-04  -S.2738E-04  0.0000E+00  0.0008E*00 

1.41 49E-I2  1 . 8443E-I1  1.I443E-01  -7.9559E-03  -2.I542E-I3  -2.0542E-03  1.9200E-01  1.337SE-I2 

I.NKErH  -1.7467E-81  1.7467E-01  8.0888E*00  3.4431E-03  -3.4431E-03  0.0000E+00  0.0800E*00 


•  *  O  O  \  v  s  - 1 


r 


Appendix  B4 


nodal  Decomposition  lor  Model  S20R5 


Eigenvalues 


-1.5186E-02  -3.9975E-82  -5.7097E-02  -1.M87E-I1  -1.3163E-01  -1.31A3E-01  -3.2738E-01  -3.523AE-01 
0.0000E+00  0.0800E+00  0.0000E+00  0.000BE+00  3.30B6E-01  -3.3086E-01  0.0000E+00  0.8000E*00 


Eigenvectors 


2.7516E-02 

9.6999E-01 

-3. 5826E-01 

-6.944BE-01 

8.0000E+00 

0.0000E+00 

0.0000E+00 

0. 0000E+00 

-8.4165E-B2 

1.5467E-01 

-2. 9S29E-81 

-6.6735E-01 

B. 00888+08 

Ba  00808  +80 

0.080BE+0B 

B.00B0E+BB 

9.3412E-01 

-1.B540E-01 

7.7235E-01 

-2.6461E-0I 

B. 0B0BE+B0 

0.00008+00 

0.8000E+00 

0.0000E+00 

-3.6414E-03 

4.9216E-04 

-5.2264E-04 

-3.5415E-03 

8.8000E+80 

0.0000E+00 

0. 00008+00 

0.0080E+0B 

5.8466E-03 

-9.38S7E-04 

2.9496E-03 

-3.2382E-04 

8.8800E+00 

8.BB80E+B0 

0.0000E+00 

0.0000E+00 

9.0608E-04 

-1.0715E-03 

1.73S0E-03 

3.7323E-03 

0.0B00E+00 

0.B008E+0B 

0.0000E+00 

0.0000E+00 

-2.9824E-02 

-2.5B64E-03 

-6.4512E-03 

4.6243E-02 

B.B004E+BB 

fl.BBBBE+BB 

fl.BBBBE+BB 

B.BBBBE+BB 

-3.4436E-01 

2.8443E-02 

-5.2797E-02 

-1.2540E-02 

0. 0000E+00 

0. 00NE+0B 

B.0000E+00 

0.B000E+00 

2.3784E-02  2.3784E-02  2.257AE-02  -1.3B93E-01 
1.7398E-02  -1.7398E-02  0.0800E+00  B.0000E+00 

2.3921E-01  2.3921E-01  2.9231E-01  9.6823E-01 
1.3985E-B1  -1.39B5E-B1  0.B008E+0B  I.008BE+88 

B.4294E-B1  B.4294E-01  -9.3751E-01  -1.2439E-B1 
-3.4306E-01  3. 430&E-B1  B.0000E+0B  0.0000E+00 

-1.B90AE-B2  -1.B906E-B2  3.7359E-B2  3.5148E-02 
-1.B192E-BI  1.B192E-B1  l.0888E*88  0.0000E+00 

5.5B02E-B3  5.5002E-B3  4.3709E-02  1.62208-02 
3. 7303E-04  -3. 7303E-04  0.00008+00  0.00008+00 

1.5456E-B3  1.5456E-03  3.3164E-03  1.9053E-B2 
1.0937E-03  -1.B937E-B3  0.0000E+00  0.00808+00 

-2. 5487E-B1  -2.54B7E-B1  -1.2756E-01  -8.4999E-02 
1 . 3604E-B1  -1.3604E-01  B.B000E+B0  0.0000E+00 

9.7445E-83  9.7445E-B3  -1.2465E-01  -4.7581E-02 
*3. 2606E-03  3.2606E-03  0.00008+01  0.0000E+00 
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Appendix  B4  (cent  d) 


Model  S20R5 


Control  lability  Matrix 

4.8721E-01  -4. 4B81E-S3  -1.3048E-02 
1.2372E-B9  -5.S447E-18  -4.4675E-18 

-9.4814E-01  -2.1125E-01  -3.4301E+00 
-5.3563E-10  5.5416E-10  -S.841SE-18 

-2.7475Et00  -5.9278E-01  -1.1337E+00 
-8.4954E-1I  1.5149E-09  5.2273E-10 

8.8974E-01  2.5212E-01  -3.4319E+00 
-3.4SS5E-18  -4.6942E-10  -1.0160E-09 

1.3e05E-01  1.8850E-03  -2.7789E-01 
1.2245E-01  -2.6434E-02  1.046BE-01 

1.36B5E-01  1.8858E-83  -2.7789E-01 
-1.2245E-I1  2. 6434E-02  -1.0460E-01 

-4.3970E-I1  2.0629E-0I  4.0265E-01 
-1. 1865E-11  -3.37I7E-U  1.2161E-10 

7.4635E-02  -4.6328E-02  -9.4798E-01 
7. 1173E-10  3.0513E-I2  -4.8438E-10 


Observability  Matrix 

9.3412E-I2  -1.8540E-02  7.7235E-02  -2.6461E-02  8.4294E-02  8.4294E-02  -9.3751E-82  -1.2439E-02 

0.0800E*0B  I.HHEvlI  l.llliE  M  I.HHE4I  -3.4306E-I2  3.4304E-I2  I.008BE+M  l.UHE^i 

9.0608E-04  -1.0715E-03  1.7350E-03  3.7323E-03  1.5456E-03  1.5454E-I3  3.3164E-03  1.90S3E-02 

I.M0E»N  I.IINE^IB  I.NIIE»II  0.0000E*0I  1.0937E-03  -1.0937E-03  0.0000E+00  I.NiiE*U 

3.4434E-01  2.8443E-I2  -5.2797E-I2  -1.254IE-I2  9.7445E-03  9.7445E-03  -1.2465E-01  -4.7581E-I2 

0.0M0E*00  I.W09E+M  0.I000E+00  0.III0E+00  -3.2406E-I3  3.2A06E-I3  I.0000E+0I  0.0000E*00 


Appendix  B5 


v 

[V 


-1.3780E-B2 

0.0008E+00 


-5.2016E-02 

0.0000E+00 


2.360SE-01  -9.6049E-01 

0.0000E*0i  0.0B00E*80 


LA 

K 


-2.2BS5E-01  -2.1244E-01 
8.8000E*00  0.0000E*00 


3.91B9E-01 

0.  MB9E+3# 


1.7733E-01 

B.B808E+B0 


-1.2734E-B2  -8.4670E-04 
i.0000E+08  0.8000E+00 


4.7857E-03 

0. 0000E*00 


8.7070E-04 

0. 0000E>00 


1.8472E-03  1.409BE-03 

8.0080E+00  8.B8B0E+00 

-2.9714E-02  9.8925E-03 

0.B000E+8I  0. 8000E>88 


nodal  Decoapositian  for  Model  525R5 


Eigenvalues 


-6.5581E-02  -1.3479E-01  -1.5137E-BI  -1.5137E-01  -4.2265E-01  -4.2265E-01 
a.0000E>88  0.0000E+00  2.8421E-81  -2.8421E-01  1.3732E-02  -1.3732E-02 


Eigenvectors 

-8.2586E-01  -6.2030E-01  6.5875E-03 

0.00B0E+00  0.0000E+00  3.9627E-02 

-3.7152E-01  -6.08 40E-01  3.1670E-02 
8.0800E+B0  0.0008Et00  2.4393E-01 

4.2192E-01  -4.B881E-01  9.3640E-01 

0.0080E+B8  0.8B00E*8B  -5.0677E-02 

-1.4691E-03  -B. 1065E-03  1.9B22E-02 
0.0000E+00  0.0000E+00  -7.B97BE-02 

1.7415E-03  -6.59B9E-04  4.2094E-03 
B.0800E+08  0.0000E+00  2. 2436E-03 


6.5875E-03  -1. 1512E-01  -1.1512E-01 
-3.9627E-02  2.e0e>2E-02  -2.6062E-02 

3. 1678E-02  6.5027E-B1  6.5027E-01 
-2.4393E-01  5.5980E-01  -5.5980E-01 

9.3648E-01  6.9B53E-02  6.9353E-02 
5.0&77E-B2  -4.6114E-B1  4.6114E-01 

1 . 9B22E-02  3.2410E-02  3.2410E-02 
7.097BE-02  5.32B1E-02  -5.3291E-02 

4.2094E-03  8.5510E-03  8.5510E-03 
-2.2436E-83  3.4444E-02  -3.4444E-02 


2.2594E-B3  3.7002E-03  1.3B49E-03  1.3B49E-B3  1.352BE-02  1.352BE-02 

0.8080E+08  0.0000E+00  1.3322E-03  -1.3322E-B3  9.9004E-03  -9.9004E-03 

1 . 6055E-B2  7.SBB5E-02  -2.23B1E-01  -2.2381E-01  -5.4217E-02  -5.4217E-02 

0.0000E+00  0.B080E+B0  S.2374E-B2  -5.2374E-B2  -1.327BE-01  1.3278E-81 


-3.8879E-01  -2.7826E-02  -4.0591E-02  -1.7234E-02  1.3502E-02  1.3502E-02  -2.0993E-82  -2.0993E-82 
8. B008E+00  8.0B00E+00  0.I8B8E+0B  S.0808E+08  -7.3380E-04  7.338BE-B4  -7.6A33E-02  7.6633E-82 


Appendix  B5  (cont'd) 


•;?1 

>& 


Model  S25R5 


Control  lability  Matrix 

6.2761E-I1  -6.4231E-B3  4.SS21E-I2 
-4.7501E-89  -3.0449E-1B  -B.7A41E-89 

4.7783E+00  1.8454E*0B  7.0500E*00 
-1.89B3E-B9  -3.94BBE-09  -2.8645E-8B 

-7.1170E+00  -1.5989E+B8  -3.B872E+00 
1.5541E-B9  5. 1546E-09  3.481BE-0B 

2.1383E+8B  4.7327E-B1  -5.0141E+00 
-1.8844E-09  -4.6893E-10  -5.S8B1E-09 

3.2985E-01  7.4203E-83  -8.3275E-0! 
2.4104E-81  -5.6660E-B2  4.1824E-B1 

3.2985E-81  7.42B3E-83  -8.8275E-01 
-2.4104E-01  5.6660E-02  -4.1B26E-01 

3.9653E-B1  -1.79B1E-01  -1.1051E+80 
5.3461E-B1  -2.1888E-81  -3.4475E-81 

3.9653E-01  -1.7981E-BI  -1.1051E+B8 
-5.3461E-01  2. 1B8BE-81  3.6675E-01 


Observability  Matrix 

B.91B9E-B2  1.7733E-B2  4.2192E-B2  -<  888 IE -82  9.344BE-B2  9.3440E-B2  4.9853E-03  4.9B53E-03 

8.8888Er08  B.80B0E+B8  B.BB08E+BB  B.0BB8E+BB  -5.0677E-03  5.8677E-83  -4.4U4E-02  4.6114E-02 

1.8472E-B3  I.4898E-03  2.2394E-83  3.7002E-03  1.3S49E-B3  1.3B49E-83  1.352BE-82  1.3520E-82 

0.8880E+B8  8. 18«BE>84  B.M88E+80  0.0080E+88  1.3322E-IJ  -1.3322E-I3  9.9004E-83  -9.9884E-03 

3.8079E-01  -2. 7824E-02  -4.I591E-02  -1.7234E-02  1.3502E-02  1.3502E-02  -2J993E-82  -2.0993E-82 

MHK*M  B.0O00E+I8  0.0808E+00  0.I00IE+88  -7.3380E-S4  7.3388E-B4  -7.4633E-02  7.6633E-82 


-  127  - 


—  w-  W.W*,  v.  r+w\ 

Ll» *  ■  *a  *  •  "  m  •  ■  V*  5  i  1  *  •  « 


MV  VWVV  VVV  v'v--.*  *  r\  • 


^  4' 


& 


•L 

vS] 


Sq 

M 


?i  ^  il 


y 

r? 


r,; 

,  i’  '■„ 

Al 


\t‘ 

v  ■■» 


Aopendix  B6 


Hulti variable  Zeros  for  all  Model s 


Node]  S5R5 

1.5030E+07  1.0854E+07  1.5425E+03  -8.4543E-03  -4.2165E-02  -4.2165E-I2  -5.3654E-02  -1.5421E+I3 

I.IINMi  0.0909E+09  0. 0000E+00  0, 0B09E-T00  3.9774E-0I  -3.9774E-01  9. 0090E>00  0.0080E+09 


Model  S10R5 

9.0004E+10  7.0271E+04  -1.4778E-02  -7.4204E-02  -7.4204E-02  -9.6508E-02  -3.4771E+04  -3.4771E*04 
0.0000E+00  0.0000E+90  0.0000E+09  3.8722E-01  -3.8722E-91  0.0000E*09  6.0715E+04  -6.0715E+04 


Node]  S20R5 

5. 1882E+06  2.5141E+04  -2.8452E-02  -1.4188E-01  -1.4188E-01  -2.0274E-01  -2.5037E+04  -2.0296E*09 

9. 0000E+00  0.9009Et00  0, 9000E+90  3.2942E-01  -3.2942E-01  0.0000E>00  0.0000E+09  0.0000Et00 


Model  S25R5 

4. 1257E+02  -3.5002E-02  -1.7204E-01  -1.7204E-01  -2.5599E-01  -1.0904E+03  -1.0904E+03  -1.2894E+04 
-1.2804E+03  9. 9099E+90  2.6989E-01  -2.6989E-01  0.9000E*00  7.9931E+02  -7.9931E+92  0.0000E+00 


p  "Me 


Appendix  C 


Sain  Matricaa  and  Propartiaa  of  tha  Cloaad  Loop 

Plant 
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Appendix  Cl 


5.0993E-01 
5. &308E-01  - 

1. 1576E-02 
S. 6803E-01 

4. 1 466E-03  - 
5. 3033E-02  - 


8. 1776E-01 
1.3793E+00  - 

8. 5743E-02 
1 . 6640E+00 

1 . 8005E-02  - 
1 . 3045E-01  - 


1 . 4466E+00 
1 . 7005E+00  - 

3. 0242E-01 
2. 3851E+00 

7. 0806E-02  - 
5.2001E-02  - 


1.71 70E+00 
1 . 7155E+00  - 

4. 1451E-01 
2. 3852E+00 

1 . 0841E-01  - 
3. 4931E-03  - 


Kalman  Filter  Sain  Matrices 


Model  SSR5 


1. 1576E-02 

4. 1466E-03 

1 . 9263E-03 

-6. 9896E-03  -5.6557E-01 

4. 7120E+01 

-2.3271E-01  -7. &189E-01 

-7 . 7504E+00 

3. 271 7E-01 

-1.5101E-03  -4. 5847E-04 

- 1 . B951E-02  -b. 9427E-01 

2. 6916E+01 

4.2244E-01  -4.4450E-01 

5. 2359E+00 

1 . 5101 E— 03 

1 . 5507E-01 

1 . 9535E-03 

1.3183E-02  -1.1 124E-02 

5. 8194E-01 

-4.7611E+00  -2. 9808E-02 

-8. 5289E-02 

Model  S10R5 

8. 5743E-02 

1.8005E-02 

1. 7876E-03 

-8. 3415E-03  -6. 2216E-01 

3. 0833E+01 

-5. 9833E-01  -7. 1770E-01 

-7. 4749E+00 

5. 3248E-01 

-5. 2855E-03  -5.8788E-04 

-2.11 42E-02  -b. 5491E-01 

1.7290E+01 

8.3463E-01  -3. 7143E-01 

6. 1774E+00 

5.2855E-03 

2.831 9E-01 

2.0265E-03 

1 . 3008E-02  -1 . 5523E-02 

1 . 0776E+00 

-4. 8471E+00  -2.839BE-02 

-2. 5462E-01 

Model  S20R5 

3. 0242E-01 

7.0806E-02 

2. 0844E-03 

-9.0445E-03  -6. 6804E-01 

1 . 9263E+01 

-1 . 3559E+00  -6. 3832E-01 

-b. 9217E+00 

9.2461E-01 

-4. 6849E-02  -3.7150E-04 

-3.2031E-02  -6. 0270E-01 

8. 5952E+00 

1 . 6947E+00 

1 . 7 106E-01 

6.9342E+00 

4. 6849E-02 

5. 6766E-01 

2. 181 0E -03 

1 . 4902E-02  -7.3951E-03 

2. 4884E+00 

-5. 0834E+00  -b.  6319E-02 

-1 . 2386E+00 

Model  S25R5 

4. 1451 E— 0 1 

1 . 0841E-01 

2. 4834E-03 

- 1 . 2226E-02  -6. 782BE-01 

1 . 6714E+01 

-1 . 7102E+00  -6. 21 42E-01 

-b, 6887E+00 

1 . 0844E+00 

-9. 1668E-02 

1.11 44E-04 

-3. 9246E-02  -5.8804E-01 

6. 7242E+00 

2. 1371E+00 

3. 6577E-01 

7. 01 79E+00 

9. 1668E-02 

7, 0799E-01 

2.2573E-03 

1 . 6466E-02  -1. 1536E-03 

3. 3726E+00 

-3. 3772E+00  -1 . 3167E-01 

-2. 0380E+00 

-r«T»n 


Appendix  C2 


Control  Gain  Matrices 


Model  S5R5 


-1 . 0018E-02 
- 1 . 47B6E-01 
-1 . 3888E-01 
-2.0376E-02 
-2.5016E-01 
3.6744E-01 
2 . 0&08E-04 
1. 1480E-04 
4. 3158E-03 
-4.6205E-04 
2. 3 1 12E-03 


-5. 9555E-02 
-4.9755E-02 
-9. 8459E-01 
2. 1 793E+00 
-2.2573E+00 
4.3158E-02 
1.2104E-03 
1 . 2507E-04 
3. 31 25E-02 
1 . 552SE-02 
8. 1 370E-03 


-5.0536E-02 
4. 0136E-02 
-2. 4462E-01 
3. 7016E-02 
-6. 4614E-01 
2.31 12E-02 
b. 4495E-04 
1 . 2895E-04 
8. 1370E-03 
4. 8522E-03 
4 . 0875E-02 


-6. 3100E-03 
-1. 170AE-01 
-1. 1104E-01 
-4. 2047E-02 
-4 . 0679E-0 1 
4. 7272E-01 
9.9092E-04 
3.8495E-05 
6.5458E-03 
-3.2026E-04 
2. 1223E-03 


Model  S10R5 

-2. 9297E-02 
-1.0929E-01 
-6. 8330E-01 
3. 6708E+00 
-3. 2997E+00 
6.5458E-02 
5.0252E-03 
-9. 8181E-04 
4. 2779E-02 
4. 5913E-02 
1 . 0M8E-02 


-2. 1827E-02 
-2.3457E-03 
-1.6012E-01 
4. J904E-01 
-9. 9766E-01 
2. 1223E-02 
1.8244E-03 
-1.6970E-04 
1. 0618E-02 
1 . 5766E-02 
S. 1 121E-02 
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Appendix  C2  (cont'd) 


Control  Gain  Matrices 


-3. 0768E-04 
-9 . 957BE-02 
-5. 2549E-02 
-4.0965E-01 
-4. 739 1 E— 02 
b.  U00E-01 
3. 8992E-03 
-5. 6741E-05 
5 . 4639E-03 
-3 . 9308E-02 
8. 5040E-04 


Model  S20R5 

1 . 8691E-02 
-3. 1810E-01 
-4. 0681E-01 
6. 7546E+00 
-4. 6661E+00 
5. 4639E-02 
2. 8093E-02 
- 1 . 0044E-02 
5. 0574E-02 
1 . 6695E-01 
8. 1331E-03 


-2. 6678E-03 
-4.3207E-02 
-5. 9456E-02 
1 . 5694E+00 
-8. 8520E-01 
8 . 5040E-03 
3. 5058E-03 
- 1 . 5228E-03 
8. 1331E-03 
3. 6961E-02 
5. 8230E-02 


2.3797E-03 
-9. 5498E-02 
-3. 4741E-02 
-1. 1024E+00 
3. 4893E-01 
7. 7970E-01 
3. 2007E-03 
9. 6791E-04 
3. 5628E-03 
-1. 1216E-01 
4. 7694E-03 


Model  S25R5 

3.3659E-02 
~4. 4344E-01 
-3, 1 467E-01 
8. 0306E+00 
-4. 6157E+00 
3. 5628E-02 
4.8469E-02 
-1 . 7839E-02 
4. 6351E-02 
2 . 4663E-0 1 
2.2212E-04 


-3.7276E-03 
-5.8065E-03 
-5. 7352E-04 
1 . 5822E+00 
-4. 8207E-02 
4.7694E-02 
-4.5751E-03 
1 . 5633E-04 
2. 2212E-04 
1.7557E-02 
5. 7978E-02 


Poles  o f  the  Closed  Loop  Model 


-7.6439E-01  - 
-3.9014E-01  - 
-3. 5318E-01  - 


7.7A16E-01  - 
■3.4961E-01 
0. 0000E+00 


3. 3022E+02 
-2. 8652E-02  - 
-1.4188E-01  - 


6. 5409E+02  - 
0. 0000E+00 
3. 2942E-01 


Real  Part 


7. 6489E-01 
1.3606E-01 
2. 0624E-01 
2. 8649E-02 


-6. 2486E-01 
-1 . 3606E-01 
-l . 2874E-01 
-1 . 303 1 E -02 


-3. 1063E-01 
-1. 5198E-01 
-5.7009E-02 
-1 . 3031E-02 


-3. 1063E-01 
-1 . 5198E-01 
-5. 7009E-02 
-4. 6457E-02 


leapinary  Part 


7. 7616E-01 
3. 3397E-0I 

0. 0000E+00 
0 . 0000E+00 


0. 0000E+00 
-3. 3397E-01 
0. 0000E+00 
b. 1 984E-03 


4. 5795E-01 
3. 2898E-01 
5. 1 976E-02 
-b. 1984E-03 


-4. 5795E-01 
-3. 2898E-01 
-5. 1976E-02 
0. 0000E+00 


Zeros  o i  the  Closed  Loop  Model 


3. 3022E+02 
9.091  IE-02 
2.0274E-01 
5. 3055E+01 


Real 

1 . 7580E+02 
-9 . 09 1 1 E — 02 
-3. 2526E-01 
-8.8860E+01 


Part 

5. 2209E+01 
-1 . 3101E-01 
-3. 4898E-01 
-8. 8860E+01 


-1.7678E-02 
-1 . 3101E-01- 
-4 . 2233E-0 1 
-8. 1 2&BE+02 


Imaginary  Part 


6. 5409E+02 
2. 2508E-02 

0. 0000E+00 
0.  0000E+00 


0. 0000E+00 
-2. 2508E-02 
0.0000E+00 
1.5280E+02 


0. 0000E+00 
3. 0736E-01 
0. 0000E+00 
- 1 . 5280E+02 


3. 0012E-03 
-3 . 0736E-01 
5. 2839E+01 

0. 0000E+00 


-3.9014E-01 
-3. 2827E-01 
-5.7206E-02 


3. 4961E-0 1 
0. 0000E+00 
0.0000E+00 


-1.7678E-02 
1 . 4188E-01 
-4.2234E-01 


-3.0012E-03 
3. 2942E-01 
-5.2639E+01 
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Appendix  D 


Gain  Scheduling  Algorithm  and  Coefficient a 


The  program  in  appendix  D1  is  written  in  Fortran  77  for 
the  Microsoft  compiler.  It  should  run  on  any  computer 
having  a  Fortran  77  compiler  with  modifications  required 
only  for  the  input /output  statements. 
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Aooendix  D1 


Fortran  Program  to  Produce  the  Gain  Scheduling  Coefficients 


program  matlsq 

parabolic  least  squares  fit  with  order  of  polynomial  taken 
as  input  from  the  console  (max  order  -  3) 
adapted  froe  the  Fortran  prograe  by  Alan  R.  Miller 
this  version  handles  matrices  up  to  11  x  11 

integer  maxr,maxc,lines,nrow,ncol  .matr, mate, max  sir 
real  x (5) ,y«at < 1 1 , 1 1 , 5) , y (51 
real  coef ( 1 1 , 1 1 , 5 ) . correl  , c  (4 ) 
character*l  answer 

character*15  fname(5),  dskfile.  prnfile 

character*60  title 

common  /sizes/  nr ow , ncol , matr . mate 

common  /files/  fname 

data  maxr .rnaxc .maxsi z /5»  4,  11/ 

wr i te (* . ' ( A) ' )  '  Specify  title  for  this  run:  (60  chars)  ' 
read(*,103)  title 

call  input (x,ymat.maxr,maxc,aaxsiz) 

Now  do  the  least  squares  fit,  storing  the  ncol  coeff  s  for  the 
i,j  term  in  the  matrix  coef f ( i , j , k )  for  k=l  tD  ncol.  Store  the 
correlation  coefficient  in  coef f ( i , j , ncol  + 1 ) . 

do  20  isl,matr 
do  20  j=l,matc 
do  10  k  =  l ,nrow 

v <k)  =  ymat(i , j.k) 
call  linfit<x,y,c, correl, maxr.maxc) 
do  IS  k=t,ncol 

coef (i,j,k)  =  c(k) 
coef (i , j ,ncol*l )  »  correl 
continue 


Appendix  D1  (cont 'd) 


Fortran  Program  to  Produce  the  Gain  Scheduling  Coefficients 


c 

c  Determine  what  type  of  output  is  desired 
c 

write!*, '  ( A \ ) ' )  '  Write  coefficients  to  output  device?  (Y/N) 
read  <*, 100) answer 

if(answer  .eq.  'y‘  .or.  answer  .eq.  'Y')  then 

write!*, ' (A\) ' )  '  Disk  file  for  Draper?  (Y/Ni  ' 
read (*, 100)answer 

if(answer  .eq.  'Y'  .or.  answer  .eq.  'v')then 

write!*, ' (A)  ' )  '  Enter  the  output  filename:  ' 
wr i te ( * ,  ’  ( A)  ' )  ’  (example:  ' ' b : coef f . dat  '  ' )  ' 
read (*,*)  dskf i le 

wr i te (* ,  '  ( A)  ' )  '  Writing  output  file. . .' 

open  (7 ,  FILE  =  dskf i le,  STATUS  =  'NEW') 
wr i te < 7, 104)  title 

wr i te ( 7, ' ( A)  ' )  '  (first  coeff  is  the  constant  term) 
do  30  i*l,matr 

do  30  jsl.matc 

30  write(7, 105) i , j , (coef ( i , j ,  k ) ,  k  =  l,  ncol) 

write(*. ' (A)  ' )  '  ' 
close(7) 
endi  f 

write!*, ' < A \ ) ' )  '  Disk  file  for  printing?  !Y/N)  ' 
read (*. 100)answer 

iftanswer  .eq.  'Y‘  .or.  answer  .ea,  *v')then 
write!*, ' (A) ' >  '  Enter  the  output  filename: 
write!*,  '  (A)  ' )  ’  (example:  ' 'dicoef f . 1st ' ' )  ' 
read<*,*)  prnfile 

wri te (*,  '  (A)  ' )  '  Writinq  output  file.. . ' 

open (8,  FILE  *  prnfile,  STATUS  =  NEW) 
write(8,104)  title 
write(8, 101) 
do  35  i=l,matr 

do  35  j=l,matc 

35  write(8,102)  i , j,coef (i , j,ncol+l) , 

1  (coef ( i . j . k ), k=l , ncol ) 

cl  ose (8) 
endif 
endif 


Appendix  Di  (cont  d) 


Fortran  Program  to  Produce  the  Gain  Scheduling  Coefficients 


100  format(al) 

101  format!'  eleaent ' ,3x , 'correlation  '  ,3x , 

1  'coefficients  flow  order  to  high  order) '/) 

102  foreatf'  ( ' , i 2, ' , ' , i2, ' ) ' , 3x , el0. 4 , 3x ,5 ( e 1 0 .  4 ,2x ) ) 

103  format!a60) 

104  format ( '  ' ,a60) 

105  fornat(2x,i2,lx,i2,lx,el2.6,2x,el2.6,2x,el2.6,2x,el2.6> 
end 

c 

subroutine  input(x,ymat,maxr,maxc, maxsiz) 
integer  nrow, i ,maxr , aaxc ,ncol ,eatr , mate , max  si z 
real  x (5) ,  vmat (11,11,5) 
character*15  fname(5) 
common  /sizes/  nr ow , ncol , eatr , mate 
common  /files/  fname 
c 

5  wri te (*, ’ ( A \ ) ' )  '  Order  of  polynomial  to  use?  ' 

read (*,*)ncol 

iffncol  .gt.  maxc-l)qoto  5 
iffncol  .It.  1)  stop 
ncol  *  ncol+1 

10  write!*, ' (A\) ' )  '  Number  of  speeds?  ' 

read (* , *) nr ow 

iffnrow  .It.  ncol  .or.  nrow  .gt.  maxr)goto  10 
write!*,  '  (A) ')  '  Enter  the  soeeds  and  correspondi np  filenames: 
write!*, ' (A) ')  '  (example:  20.0  '  'b: abohe. s20 ' ' ) ' 

read!*,*)  (x (i ) .fname d )  ,  i=l,nrow> 

15  wr i te (*, ' ( A \ ) ' >  '  Specify  matrix  dimensions:  (row, col)  ' 

read!*,*)  matr.matc 

iffmatr  .gt.  maxsiz  .or.  mate  .gt.  maxsiz)  goto  15 

write!*,  '  (A)  ' )  '  Reading  input  matrices . ' 

do  20  i=l,nrow 

openfi,  FILE  *  fname(i)) 

read(i,100)  ( ( j , k , ymat ( j ,  k ,  i > ,  m=l,matc),  n=l,matr) 
closed ) 

20  continue 

write!*, 101) 
return 

100  f or mat (Ix,i2,lx,i2,lx,el6.10,lx,i2,lx,i2,lx,el6.10, 

1  lx , i2, lx , i 2, lx ,e!6. 10) 

101  format('0'> 
end 
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Appendix  Dl  (cont'd) 


Fortran  Program  to  Produce  the  Gain  Scheduling  CoeFFicients 


subroutine  linFit (x , y , coeF , cor ,»axr , eaxc) 
logical  error 

integer  nroM,ncol,i,j,eaxr, eaxc ,eatr , eatc 
integer  index (5,5) ,nvec 
real  x (1 ) , y ( 1 ) ,coeF  ( 1 ) 
real  a<5,5) ,x»atr (5,5) 
real  suey,su«y2,xi,yi,yc,res,cor,srs 
coaaon  /sizes/  nrow,ncol ,*atr ,aatc 
data  nvec/1/ 
c 

do  10  i =1 , nrow 
x  i  =  x  ( i ) 
xoiatr  < i , 1 ) =1 . 0 
do  10  j=2.ncol 

x*atr(i,j)=  xmatr  (i ,  j-1 ) *xi 
10  continue 

call  squarelxeatr.y.a.coeF.maxr.maxc) 
call  g aussjla, coeF, index, nvec. err or, maxc) 

SUfflV  =  0.0 
suffly2  =  0.0 
srs  -  0.0 
do  20  i*l.nroN 
yi  =  y  (i ) 
yc  =  0.0 

do  15  j  =  1,  ncol 

15  yc  *  yc  ♦  coeF(j)  *  xeatr(i.j) 

res  *  yc  -  yi 

srs  =  srs  +  res*res 

su*v  =  suey  ♦  yi 
su»y2  *  su«v2  ♦  yi  *  yi 
20  continue 

c  take  care  oF  the  case  where  the  v(i)  are  all  zero 

c  this  is  guaranteed  to  produce  zero  coeFFicients  so  cor  *  1.0 
iFisuay  .eq.  0.0)  then 
cor  *  1.0 

else 

cor  -  sqrt  (1 . 0-srs/ (suey2-su»y*su»y/nro*i) ) 
endif 
return 
end 


7 


Appendix  Dl  (cont'd) 

Fortran  Program  to  Produce  the  Gain  Scheduling  Coefficients 


subroutine  square (x , y , a ,g ,eaxr , maxc > 
integer  nr o*,ncol,i,k,l,matr, mate, maxr, maxc 
real  x  (maxr , eaxc ) ,y(maxr),a(maxc, maxc), g( maxc) 
common  /sizes/  nrow,ncol ,*atr ,matc 


do  40  k»l,ncol 
do  20  1=1, k 
a  (k ,  1 ) =0.0 
do  10  i=l,nrow 

a (k ,1 )=a(k , 1 ) +x  <i , 1 )*x (i ,k> 
i f ( k  . n  e .  1)  a  ( 1 , k )  a  a  ( k  .  1 ) 
continue 
continue 
g(k)  =  0.0 
do  30  i=l,nroM 

g  ( k  1  *g  (k )  +y  ( i  1  #x  (i  •  k  1 
continue 
continue 
return 
end 


subroutine  oaussj(b.w, index . nvec , error . naxc ) 
logical  error 

integer  nrow, i , j , k , 1 ,nvec ,matr .mate .naxsi z .eaxc 
integer  i row, i col , 1 1 , i ndex (eaxc. 3) 
real  b(maxc, l),**(aaxc,l), big, sum, t, pivot, determ 
common  /sizes/  nrow.ncol ,matr .mate 

error  =  .false, 
n  *  ncol 
do  10  i =1 ,n 
index (i ,3) *0 
continue 
determ  =  1.0 
do  90  i*l,n 
big  *  0.0 
do  20  j ■ l . n 

if  (index (j, 3)  .eq.  llgoto  20 
do  IS  k«l,n 

if  (index (k, 3)  .gt.  llgoto  199 
i f (index (k ,31  .eq.  1)  goto  15 
if <abs<b(j,kJ J  .le.  biglgoto  IS 
iron  ■  j 
icol  ■  k 

big  •  mbt(b(j,k) ) 
continue 
continue 
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Appendix  D1  (cont 'd ) 

Fortran  Program  to  Produce  the  Gain  Scheduling  Coefficients 


index (icol ,3)  =  index (icol ,3)  ♦  1 

index (i ,2)  *  irow 

index (i , 1 )  *  icol 

if(irow  .eq.  icollgoto  40 

determ  =  -deter* 

do  25  1  *1 ,  n 

call  swap (b (i row, 1 ) ,b (i col , 1 ) > 
continue 

i  f (nvec  . eq.  0) qoto  40 
do  30  1*1, nvec 

call  swap (w ( i row, 1 ), w ( i col  ,  1 ) ) 
continue 

pivot  *  b ( icol , i col ) 
determ  *  determ*pivot 
b  ( l col , i col ) *  1.0 
do  45  1=1 ,n 

b(icol,l)=b(icol,l)/pivot 

continue 

i f  (nvec  . eq.  0) goto  60 
do  50  1*1, nvec 

w(icol,l)  =  w(icol,l)/pivot 
continue 
do  80  1 1=1 , n 

if  (1 1  .eq.  icollgoto  80 
t  =  b  < 1 1 , icol ) 
bill, icol)  =  0.0 
do  65  1=1, n 

b  ( 1 1 , 1 )  =  b  (1 1 , 1 1 -b (icol ,  1 )  *t 
continue 

i f  (nvec  . eq.  0) goto  80 
do  70  1=1, nvec 

w(ll.l)  =  w ( 1 1 , 1 ) -w ( i col ,  1 1 *t 
continue 
continue 
continue 
do  120  i=l,n 
1  *  n  -  i  4-1 

i f  (index (1 , 1 )  .eq.  index  (l  ,2) Igoto  120 
irow  ■  index (1,1) 
icol  *  index (1 .2) 
do  110  k  *  1 ,n 

call  twap(b(k, irow) ,b(k, icol 1 ) 
continue 
continue 
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Appendix  Di  (cont'd) 


Fortran  Progra*  to  Produce  the  Gain  Scheduling  Coefficients 


do  130  k  =  l,n 

if (index (k, 3)  .ne.  llgoto  199 
130  continue 

return 

199  write!*, 999) 

error  =  .true, 
return 

999  foreat!'  ERROR  -  natrix  singular') 

end 
c 

subroutine  swap(a.b) 

real  a, b, hold 

hold  =  a 

a  *  b 

b  =  hold 

return 

end 


Appendix  D2 


A-BG-HC  matrix  coefficients 

coefficients 


element 

correlation 

(low  order  to  high 

order ) 

( 

1, 

1) 

. 9999E+00 

1694E+00 

-. 6835E-01 

. 2494E-03 

( 

1, 

2) 

. 9996E+00 

. 6212E— 01 

1287E-01 

-.252BE-03 

( 

1, 

3) 

. 1000E+01 

. 2333E-02 

-.4777E-03 

1583E-03 

( 

1, 

4) 

. 9998E+00 

— . 2230E— 02 

. 7932E-04 

-.3585E-05 

( 

1 , 

5) 

. 9479E+00 

. 7586E-02 

1184E-03 

. 1153E-04 

< 

1, 

6) 

. 9972E+00 

. 5033E+00 

. 1461E-01 

3083E— 03 

( 

1 , 

7) 

. 9850E+00 

. 3162E+00 

2128E+00 

. 5342E-02 

( 

1, 

8) 

. 9954E+00 

. 6508E+02 

4195E+01 

. 9159E-01 

( 

1 , 

9) 

. 9999E+00 

5577E-01 

. 6878E— 01 

. 2617E-04 

< 

1, 

10) 

. 9987E+00 

. 8225E+00 

1267E-01 

. 1820E-03 

( 

1, 

11) 

. 9999E+00 

. B144E+01 

— . 7023E— 01 

. 4817E-03 

( 

2, 

1) 

. 9996E+00 

. 6212E-01 

1287E-01 

2528E-03 

( 

2, 

2) 

. 9999E+00 

9588E-01 

— . 4723E— 01 

. 3033E-03 

( 

2, 

3) 

. 9995E+00 

. 1423E-01 

-.3771E-02 

. 2736E-03 

( 

2, 

4) 

. 9966E+00 

. 7009E— 04 

. 9534E-04 

4082E-05 

< 

2f 

5) 

. 9954E+00 

. 2375E-02 

. 1945E-02 

1984E-04 

( 

2, 

6) 

. 9999E+00 

. 7592E+00 

1101E-01 

. 1B25E-03 

< 

2, 

7) 

. 9968E+00 

. 7432E+00 

3063E+00 

. 7306E-02 

< 

2, 

8) 

. 9985E+00 

3801 £+02 

.  2533E+01 

-.5170E-01 

< 

2, 

9) 

. 1000E+01 

1555E-01 

7298E-01 

. 2372E-03 

( 

2, 

10) 

. 9921E+00 

. 5869E+00 

1893E-01 

—.81 10E— 03 

< 

2* 

11) 

. 9976E+00 

4239E+01 

— .  2394E+00 

.5186E-02 

( 

3, 

1) 

. 1000E+01 

. 2333E-02 

— . 4777E— 03 

-. 1583E-03 

< 

3, 

2) 

. 9995E+00 

. 1423E-01 

— . 3771E-02 

. 2736E-03 

( 

3, 

3) 

. 9999E+00 

-.2462E-01 

— . 2536E— 01 

— . 8143E-04 

( 

3, 

4) 

. 1000E+01 

1 879E-02 

1460E-04 

— . 2144E— 07 

< 

3, 

5) 

. 9984E+00 

1373E-01 

. 1785E-03 

1159E-04 

< 

3, 

6) 

.9912E+00 

.232 IE -01 

.  8328E-03 

— . 6426E— 04 

< 

3, 

7) 

. 9227E+00 

.2195E-01 

2131E— 01 

. 8262E-03 

< 

3, 

8) 

. 1000E+01 

. 2034E+00 

.6213E-01 

. 2590E-02 

( 

3, 

9) 

. 9992E+00 

.611 1E+01 

— . 8634E— 01 

. 2780E— 02 

( 

3, 

10) 

. 9952E+00 

. 6388E-01 

— . 8520E— 02 

. 445 IE— 03 

C 

3, 

11) 

. 9999E+00 

. 4404E+00 

— .  4512E-01 

. 4371E— 02 

< 

4, 

1) 

.9998E+00 

. 1 433E-03 

8452E— 04 

-. 2997E-04 

< 

4, 

2) 

. 9998E+00 

. 4106E— 03 

2215E-03 

7572E— 04 

c 

4, 

3) 

. 9998E+00 

. 1 652E-02 

- . 8909E-03 

- . 3046E— 03 

( 

4, 

4) 

. 9998E+00 

1616E-02 

-. 1360E-02 

. 4314E-06 

< 

4, 

5) 

. 9994E+00 

-.B778E-03 

2140E— 02 

. 9642E-05 

( 

4, 

6) 

. 9907E+00 

. 3420E-01 

-.7847E-02 

. 4101E-03 

c 

4, 

7) 

. 9994E+00 

-. 7683E-03 

1833E-02 

. 778 IE— 05 

( 

4, 

8) 

. 995BE+00 

-. 4556E+00 

. 1650E+00 

-. 8371E-02 

< 

4, 

9> 

. 9967E+00 

— . 3616E+00 

5120E-01 

.61 85E-03 

< 

4, 

10) 

. 1000E+01 

•0000E+00 

. 0000E+00 

. 0000E+00 

< 

4, 

11) 

. 984BE+00 

. 7706E+00 

1587E+00 

. 2462E— 02 
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Appendix  D2  (cont'd) 


A-BG-HC  matrix  coefficients 


coefficients 


element 

correlation 

(low  order  to  high 

order ) 

( 

5,  1) 

. 1000E+01 

. 0000E+00 

. 0000E+00 

. 0000E+00 

( 

5,  2) 

. 1000E+01 

.  0000E+00 

. 0000E+00 

. 0000E+00 

( 

5,  3) 

. 9998E+00 

2143E-01 

. 8462E-02 

. 2433E-02 

< 

5,  4) 

. 1 000E+0 1 

. 3066E-03 

. 2370E-03 

. 1173E-05 

5,  5) 

. 9997E+00 

— . 5286E-02 

— . 6884E-02 

. 7536E-05 

5,  6) 

. 9998E+00 

1165E-01 

. 1008E-01 

-. 4830E-03 

5,  7) 

. 9666E+00 

. 1818E+00 

1 226E+00 

. 4071E-02 

5,  8) 

. 9735E+00 

1456E+00 

. 6308E-01 

-. 2428E-02 

5,  9) 

. 1000E+01 

. 1911E+00 

— . 5891E— 01 

-. B543E— 02 

5,10) 

. 9907E+00 

. 1214E+00 

. 25 1 BE— 02 

-. 1 306E-03 

5,11) 

. 9870E+00 

-.2566E+00 

.21 30E+00 

8000E-02 

6,  1) 

. 9998E+00 

. 1688E— 01 

-.6976E-02 

-.2069E-02 

6,  2) 

. 9998E+00 

.81 78E— 02 

— . 4408E— 02 

-. 1507E-02 

6,  3) 

. 9998E+00 

. 4192E— 07 

— . 2261E— 07 

— . 7728E-08 

6,  4) 

. 9987E+00 

— . 2842E-03 

. 1 458E-04 

-. 3137E-05 

6,  5) 

. 9952E+00 

. 4256E-02 

-. 9237E-03 

. 5646E-04 

6,  6) 

. 9998E+00 

— . 3663E-0 1 

— . 4365E— 02 

— . 6024E— 04 

6,  7) 

. 9994E+00 

— . 6944E— 01 

-. 1435E+00 

. 3339E-03 

6,  8) 

. 9997E+00 

. 3099E+00 

. 3689E+00 

— . 4837E— 03 

6,  9) 

.9990E+00 

. 6949E-01 

-. 3453E-01 

. 1872E-02 

6,10) 

. 1000E+01 

. 0000E+00 

. 0000E+00 

. 0000E+00 

6,11) 

. 8381E+00 

— . 5165E-01 

. 7623E-02 

— . 2808E-03 

7,  1) 

. 1000E+01 

. 0000E+00 

. 0000E+00 

. 0000E+00 

7,  2) 

. 1000E+01 

. 0000E+00 

. 0000E+00 

. 0000E+00 

7,  3) 

. 9998E+00 

2036E-03 

. 1237E-03 

. 4434E-04 

7,  4) 

. 9613E+00 

. 6191E— 04 

-.2197E-04 

. 8259E-06 

7,  5) 

. 1000E+01 

5701E-03 

— . 3009E-03 

-. 3688E-05 

7,  6) 

. 9782E+00 

.6719E-03 

- . 2997E— 03 

. 6828E-05 

7,  7) 

. 9998E+00 

1383E-01 

-. 1407E— 01 

-. 5896E-05 

7,  8) 

. 9847E+00 

.  7534E-01 

-.4371E-01 

. 9755E-03 

7,  9) 

. 9872E+00 

. 3083E— 0 1 

— . 1 128E-01 

. 492 IE— 03 

7,10) 

. 9907E+00 

1493E+00 

-. 3097E-02 

. 1606E-03 

7,11) 

. 9774E+00 

.  1430E-01 

- . 3452E-02 

. 1564E-03 

8,  1) 

.9998E+00 

.  4363E— 03 

1735E-03 

— . 5014E— 04 

8,  2) 

. 9998E+00 

7851E-04 

. 4234E-04 

. 1447E-04 

8,  3) 

.9998E+00 

1383E-07 

. 7456E-08 

. 2549E-08 

8,  4) 

.9998E+00 

1422E-05 

-. 1105E-05 

-. 2545E-07 

8,  5) 

.9973E+00 

2191E-04 

. 7237E-05 

-. 4019E-06 

8,  6) 

. 9964E+00 

- . 6585E— 04 

.  571 IE— 04 

9198E-06 

8,  7) 

.  9994E+00 

— . 4591E-03 

-.  U05E-02 

. 481 IE— 05 

8,  8) 

.9998E+00 

1 466E-01 

-. 1386E-01 

-. 1341E-04 

8,  9) 

.9712E+00 

— . 2680E-02 

. 3356E-03 

-. 1202E-05 

8,10) 

. 1000E+01 

.0080000 

•  WWWL  *  Ull 

e  00B0E+B0 

8,11) 

. 8166E+00 

4507E-02 

. 3768E-03 

-. 1 089E-04 
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Appendix  D2  (cont'd) 


A-BG-HC  matrix  coefficients 


element 

correlation 

(  9,  1> 

. 1000E+01 

(  9,  2) 

. 1000E+01 

(  9,  3) 

. 999BE+00 

(  9,  4) 

. 9999E+00 

(  9,  5) 

. 9999E+00 

(  9 ,  6) 

.9919E+00 

(  9,  7) 

. 9998E+00 

(  9,  8) 

. 9968E+00 

(  9,  9) 

. 9997E+00 

(  9,10) 

. 9907E+00 

(  9,11) 

. 9929E+00 

(10,  1) 

. 1000E+01 

(10,  2) 

. 1000E+01 

(10,  3) 

. 1000E+01 

(10,  4) 

. 1000E+01 

(10,  5) 

. 1000E+01 

(10,  6) 

. 9931E+00 

(10,  7) 

. 1400E+52 

(10,  8) 

-9871E+00 

(10,  9) 

. 9994E+00 

(10,10) 

. 9935E+00 

(10,11) 

. 9883E+00 

(11,  1) 

. 1000E+01 

(11,  2) 

. 1000E+01 

(11,  3) 

. 1000E+01 

(11,  4) 

. 1000E+01 

(11,  S) 

. 1000E+01 

(11,  6) 

. 8302E+00 

<11,  7) 

. 1000E+01 

(11,  8) 

. 9926E+00 

(11,  9) 

. 9999E+00 

(11,10) 

. 9996E+00 

(11,11) 

. 9979E+00 

coefficients 

(low  order  to  high  order) 


. 0000E+00 
. 0000E+00 
. 57&8E-03 
3588E-05 
-• 6080E-04 
1029E-03 
1297E-02 
. 6898E-02 
3057E-01 
. 2412E-03 
1 026E-02 
. 0000E+00 
.  0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+e0 
. 3843E-02 
. 1000E+01 
. 6741E-01 
1109E+00 
U25E-08 
. 1915E-01 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+00 
-.5870E-03 
. 0000E+00 
. 9916E+00 
. 1218E-01 
. 7924E-03 
282BE-01 


. 0000E+00 
. 0000E+00 
2201E-03 
— . S043E— 05 
— . 4190E-04 
- . 8578E— 04 
. 8453E-04 
1812E-02 
1769E-01 
. 5003E-05 
1709E-02 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+00 
-.9352E-03 
-0000E+00 
1600E-01 
. 2537E-01 
. 2765E-09 
6268E-02 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 7487E-04 
.  0000E+00 
. 2057E-02 
4383E-02 
. 1129E-02 
-.2926E-02 


. 0000E+00 
. 0000E+00 
-.6176E-04 
— . 6170E— 08 
1457E-06 
. 3225E-05 
-.2108E-04 
. 1082E-03 
. 949&E-04 
— . 2594E— 06 
-6912E-04 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 0000E+00 
. 4867E-04 
. 0000E+00 
. 8039E-03 
241SE— 02 
1371E-10 
. 1176E-03 
. 0000E+00 
. 0000E+00 
. 0000E+00 
.  0000E+00 
. 0000E+00 
2739E-05 
. 0000E+00 
1105E-03 
. S224E-03 
1651E-05 
. 6999E-04 


Appendix  02  (cont'd) 


6  matrix  coefficients 


coefficients 


el ement 

correlation 

(low  order  to  high 

order > 

( 

1. 

1) 

. 9999E+00 

. 1&94E+00 

. 6835E-01 

2494E-03 

C 

1. 

2) 

. 9996E+00 

6212E-01 

. 1287E-01 

.  2528E-03 

( 

1« 

3> 

. 1000E+01 

2333E-02 

. 4777E-03 

. 1583E-03 

( 

1. 

4) 

• 9998E+00 

. 2230E-02 

-.7932E-04 

. 3585E-05 

< 

1. 

5) 

. 9479E+00 

-.75B6E-02 

. 1184E-03 

1153E-04 

< 

1. 

6) 

. 9972E+00 

- . 50 1 9E+00 

1469E-01 

. 3090E-03 

( 

1, 

7) 

. 9850E+00 

— . 3162E+00 

. 2128E+00 

- . 5342E-02 

( 

1» 

8) 

. 9954E+00 

-.6508E+02 

. 4195E+01 

9159E— 01 

< 

1* 

9) 

. 9999E+00 

. 1518E+00 

7652E-01 

. 7580E-04 

( 

1* 

10) 

. 9987E+00 

8225E+00 

. 1267E-01 

1820E-03 

( 

it 

11) 

. 9999E+00 

- . 8058E+0 1 

. 6203E-01 

— . 2832E— 03 

< 

2, 

1) 

. 9996E+00 

6212E-01 

. 1287E-01 

. 2528E-03 

( 

2, 

2) 

. 9999E+00 

. 9588E-01 

. 4723E-01 

-.3033E-03 

( 

2, 

3) 

. 9995E+00 

1423E-01 

. 3771E— 02 

2736E— 03 

( 

2, 

4) 

. 9966E+00 

- . 7009E— 04 

- . 9534E— 04 

. 4082E-0S 

( 

2» 

5) 

. 9954E+00 

2375E-02 

1945E-02 

. 1984E-04 

( 

2, 

6) 

.  1000E+01 

-.7411E+00 

. 1023E-01 

1&4&E-03 

< 

2, 

7) 

. 9968E+00 

7432E+00 

. 3063E+00 

730&E-02 

( 

2, 

8) 

. 9985E+00 

.380 IE +02 

— • 2533E+01 

. 5170E-01 

( 

2, 

9) 

. 1000E+01 

. 1B13E-01 

. 7974E— 01 

. 2015E-03 

< 

2, 

10) 

.9921E+00 

— . 5869E+00 

. 1B93E-01 

.81 10E-03 

( 

2, 

11) 

. 9982E+00 

. 41 16E+01 

. 2580E+00 

-.5719E-02 

( 

3, 

1) 

. 1000E+01 

- . 2333E-02 

. 4777E-03 

. 1583E-03 

< 

3* 

2) 

. 9995E+00 

1423E-01 

. 3771E— 02 

- . 2736E— 03 

( 

3, 

3) 

. 9999E+00 

. 2462E-01 

. 2536E— 01 

. 8143E-04 

< 

3, 

4) 

. 1000E+01 

. 1879E-02 

. 1460E-04 

. 2144E-07 

( 

3, 

5) 

. 9984E+00 

. 1373E-01 

1785E-03 

. 1159E-04 

< 

3» 

6) 

. 9B20E+00 

-.5691E-02 

1567E-02 

. 7094E-04 

( 

3, 

7) 

. 9227E+00 

2195E-01 

. 2131E-01 

— . 8262E-03 

< 

3, 

8) 

. 1000E+01 

2034E+00 

&213E-01 

-.2590E-02 

< 

3, 

9) 

. 9955E+00 

4802E+01 

. 1219E-01 

1385E-02 

( 

3, 

10) 

. 9952E+00 

— . 6388E-01 

. 8520E-02 

— . 4451E—03 

( 

3, 

11) 

. 1000E+01 

1 199E+00 

. 2822E-01 

— . 4200E— 02 

-  143  - 


Appendix  D2  (cont'd) 


el ement 

(  1,  1) 
<  1,  2) 
(  1,  3) 
(  2,  1) 
(  2,  2) 
(  2,  3) 

<  3,  1) 

<  3,  2) 
(  3,  3) 
(  4,  1) 
(  4,  2) 
(  4,  3) 

<  5,  1) 

<  5,  2) 

<  5,  3) 
(  6,  1) 
(  6,  2) 

<  6,  3) 
(  7,  1) 
(  7,  2) 
(  7,  3) 
(  8,  1) 
(  8,  2) 

<  8,  3) 
(  9,  1) 

<  9,  2) 
(  9,  3) 
(10,  1) 
(10,  2) 
(10,  3) 
(11,  1) 
(11,  2) 
(11,  3) 


H  matrix  coefficients 


coefficients 

(low  order  to  high  order) 


correlation 

. 9999E+00 
. 9999E+00 
. 9976E+00 
.  991 1E+00 
. 9997E+00 
. 9815E+00 
. 9988E+00 
. 9975E+00 
. 9974E+00 
. 9921E+00 
. 9998E+00 
. 9848E+00 
. 9978E+00 
. 9983E+00 
. 9870E+00 
. 9903E+00 
. 9804E+00 
. 8302E+00 
. 9561E+00 
. 1000E+01 
. 8933E+00 
. 9153E+00 
. 1000E+01 
. 7781E+00 
. 9804E+00 
. 9983E+00 
. 9929E+00 
.  9931E+00 
. 9998E+00 
. 9140E+00 
.8302E+00 
. 9929E+00 
. 9979E+00 


1391E-01 
9601E-01 
8577E-01 
1806E+00 
2579E-02 
. 1236E+00 
1752E+00 
1309E+01 
- . 3206E+00 
— . 361 1E+00 
. 4895E+00 
7704E+00 
. 4290E-01 
7009E+00 
. 2508E+00 
. 3376E+00 
. 1273E-01 
. 5870E-01 
2202E-02 
. 3365E-02 
— . 7054E-02 
. 9755E-03 
- . 9334E-03 
. 2180E-02 
. 1273E-02 
. 1790E-01 
. 1014E-02 
3843E-01 
2817E-02 
— . 2804E-01 
. 5870E-02 
. 1014E-02 
. 2828E-01 


. 8202E-03 
. 7746E-02 
. 8209E-02 
. 7790E-02 
6763E— 02 
1863E-01 
. 734 IE— 02 
.7416E-01 
.  1690E— 01 
•8359E-01 
. 3388E+00 
. 15B7E+00 
7944E-01 
3402E+00 
-.2115E+00 
. 6572E-02 
. 7745E-02 
-.7487E-02 
. 4746E-03 
- . 9604E-03 
. 1677E-02 
18B0E-03 
. 4378E-03 
4215E— 03 
. 7745E-03 
. 3460E-02 
. 1712E-02 
. 9352E-02 
. 1807E-02 
. 6996E-02 
— • 74B7E— 03 
. 1712E-02 
. 2926E-02 


- . 6807E-05 
1019E-03 
1985E-03 
1782E-03 
4387E-03 
. 5326E-03 
-.6687E-04 
1 394E—02 
1708E-03 
— . 4474E— 02 
1437E-02 
-.2462E-02 
. 3686E-02 
. 7286E-02 
. 7923E-02 
. 4228E-03 
— » 2754E— 03 
. 2739E-03 
9889E— 05 
. 1104E-03 
&174E-04 
. 7340E-05 
4459E— 04 
. 1319E-04 
2754E-04 
— . 9251E— 04 
— . &928E— 04 
4867E— 03 
. 3286E-03 
-.2021E-03 
.2739E-04 
— »  6928E— 04 
6999E-04 
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